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BRDRFLW  is  the  documentation  for  a computer  program  of  the  same  name, 
\written  in  FORTRAN  and  designed  to  predict  the  behavior  of  the  surface 
stream  flowing  down  an  irrigation  border.  The  infiltration  is  assumed  a given 
function  of  infiltration  time.  Other  physical  characteristics  required  as  input 
to  the  model  are  the  border  length,  slope,  roughness,  and  downstream  boun- 
dary condition;  that  is,  whether  the  stream  is  blocked  by  an  end  check  or 
allowed  to  drain  freely  into  a drainage  ditch.  The  management  parameters— 
required  depth  of  application,  inflow  rate,  and  cutoff  time— complete  the 
physical  input  to  the  model.  Upon  entry  of  certain  numerical-solution  and 
display  parameters,  the  computer  program  solves  the  equations  governing  the 
flow  in  the  surface  stream.  The  key  solution  variables  are  the  advance,  reces- 
sion, and  runoff  as  functions  of  time,  and  the  ultimate  distribution  of  applied 
water.  These  results  can  be  displayed  in  both  tables  and  graphs.  Also  dis- 
played are  various  figures  of  merit  of  the  irrigation;  for  example,  distribution 
and  application-efficiency  parameters.  Input  can  be  made  in  either  English  or 
metric  units  or  in  dimensionless  form;  output  appears  in  all  three. 

The  governing  (nonlinear)  equations  are  comprised  of  a mass-conservation 
relation  and  either  a statement  of  equilibrium  among  pressure,  gravity,  and 
drag  forces  on  elements  of  the  surface  stream  (zero-inertia  formulation)  or 
one  of  equilibrium  between  gravity  and  drag  forces  (kinematic-wave  model, 
suitable  for  steep  slopes).  The  user  can  also  choose  a hybrid  model,  utilizing 
zero-inertia  concepts  for  advance  and  kinematic-wave  analysis  for  recession. 
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Foreword 


This  manuscript  represents  the  results  of  several  decades  of  work  on  the 
development  of  models  to  describe  the  flow  of  water  within  irrigation  borders. 
Of  the  two  basic  types  of  surface  irrigation,  borders  and  furrows,  borders  are 
the  simplest  mathematically  and  so  have  been  given  the  most  attention  ini- 
tially. The  purpose  of  surface-irrigation  models  is  to  improve  the  overall  state- 
of-the-art  of  surface  irrigation.  This  improvement  could  eventually  lead  to  the 
more  efficient  use  of  water  and  water  conservation.  The  U.S.  Water  Conserva- 
tion Laboratory  became  interested  in  surface  irrigation  modeling  in  the 
mid-1970’s  after  the  Western  Regional  Research  Project  on  surface  irrigation 
modeling  W-65  was  terminated.  This  group,  along  with  others,  had  made 
significant  advances  in  the  theoretical  development  of  border-irrigation  flow. 

A major  breakthrough  was  the  development  of  the  zero-inertia  model  by 
Theodor  Strelkoff  and  Nikolaos  D.  Katopodes  (under  cooperative  agreement 
between  the  U.S.  Water  Conservation  Laboratory  and  the  Department  of  Land, 
Air,  and  Water  Resources,  University  of  California,  Davis).  This  information 
was  published  in  1977.  The  model  had  significant  advantages  over  previous 
models  since  it  was  far  less  complicated  than  the  more  complete  hydro- 
dynamic  models  and  far  more  widely  applicable  than  the  kinematic-wave 
models. 

Versions  of  this  original  model  are  being  used  by  a number  of  researchers. 
These  versions  have  been  modified  to  handle  a variety  of  situations;  however, 
a number  of  computational  problems  have  arisen  and  the  model  as  such  has 
never  been  fully  reliable.  Accordingly,  it  was  never  fully  documented  since  an 
unfamiliar  user  may  not  be  able  to  differentiate  a good  run  from  a poor  one. 

Continual  development  has  taken  place  since  the  original  model  was  produced 
through  cooperative  efforts  with  the  University  of  California,  Davis,  and  cur- 
rently through  a contract  with  Strelkoff  (formerly  of  U.C.,  Davis).  This  ongoing 
cooperative  effort  with  Strelkoff  over  the  last  6 or  7 years  has  resulted  in  the 
mathematical  model  and  computer  program  described  herein.  This  publica- 
tion consists  of  two  theoretical  models  (the  kinematic  wave  and  a new  fully 
nonlinear  zero-inertia  or  equilibrium  model)  and  a hybrid  of  the  two.  The  user 
may  select  the  appropriate  model  for  his  or  her  particular  situation.  New  tech- 
niques have  been  used  to  perform  the  computations  and  a considerable 


amount  of  flexibility  has  been  added.  This  new  program  is  sufficiently  more 
reliable  than  its  predecessor  that  sawtooth  profiles  and  aborted  runs  can  be 
avoided  and  consistent  solutions  can  be  obtained.  This  new  program,  along 
with  its  documentation,  will  provide  a useful  tool  for  analyzing  border  irriga- 
tion flow.  It  can  be  used  to  develop  design  criteria,  to  develop  or  modify 
specific  designs,  to  evaluate  current  systems,  to  develop  optimum  manage- 
ment practices,  and  to  improve  our  understanding  of  border-irrigation  flow.  As 
such,  it  can  be  a valuable  tool  for  teaching,  for  research,  and  for  application 
to  improvements  in  border-irrigation  practices.  Those  interested  in  the  com- 
puter program  can  obtain  a copy  by  sending  a 9-track  magnetic  tape  to  me  at 
the  following  address: 


Albert  J.  Clemmens 

U.S.  Water  Conservation  Laboratory 

4331  East  Broadway 

Phoenix,  Ariz.  85040 

ATTN:  BRDRFLW 
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1.  Introduction 


1.1  BRDRFLW: 
An  Overview 


1.2  Significant  Features 


BRDRFLW  is  a computer  program  designed  to  model  the  flow  of  water  over 
the  surface  of  an  irrigation  border.  It  was  implemented  at  the  computer 
center  of  the  Lawrence  Berkeley  Laboratory,  University  of  California, 

Berkeley.  The  inflow  is  assumed  to  be  distributed  evenly  across  the  width  of 
the  border  at  its  upper  end,  and  the  transverse  slope  of  the  border  is  assumed 
to  be  zero.  Roughness  and  infiltration  characteristics  are  assumed  not  to  vary 
transversely  either,  so  that  the  resulting  flow  is  plane,  that  is,  two  dimen- 
sional. The  subject  of  analysis  is  a strip  of  field  of  unit  width  to  which  water 
is  introduced  at  a known  volumetric  rate  per  unit  width  at  the  upper  end  and 
which  possesses  known  roughness  and  infiltration  characteristics  as  well  as 
known  bottom  configuration,  length,  and  downstream-boundary  condition, 
either  a free  outflow  into  a drainage  ditch  or  ponding  behind  a dike. 

On  the  basis  of  the  physical  principles  of  conservation  of  mass  and  one  or 
another  approximation  to  the  impulse-momentum  relationship,  the  longitudinal 
variation  of  depth  and  discharge  along  the  length  of  the  surface  stream  is 
computed  at  a sequence  of  times  after  the  start  of  an  irrigation.  As  byproducts 
of  these  calculations,  the  movement  of  the  wave  front  during  stream  advance 
and  of  the  trailing  edge  of  the  stream  profile  in  recession  are  obtained.  With 
infiltration  a known  function  of  opportunity  time,  the  ultimate  post-irrigation 
longitudinal  distribution  of  infiltrated  water  depth  (volume  infiltrated  per  unit 
plan  area  of  field)  is  readily  obtained.  With  the  irrigation  requirement  known, 
the  volume  lost  to  deep  percolation  is  then  determined.  An  additional 
byproduct  of  the  profile  computations,  the  runoff  rate  from  the  end  of  the 
field,  is  integrated  over  time  to  provide  the  volume  of  runoff.  Final  computa- 
tions by  the  program  provide  the  water-application  efficiency  and  various 
distribution  parameters  stemming  from  the  inputted  conditions  of  the  irrigation. 

The  zero-inertia  border  mathematical  model  was  last  documented  in  1977  (10). 
At  that  time  a locally  linearized  set  of  algebraic  equations  was  used  to 
approximate  a mix  of  integrated  and  differential  forms  of  the  equations  of 
continuity  and  motion.  The  present  equilibrium  model  is  based  wholly  on  inte- 
grated forms.  Statements  of  mass  conservation  and  equilibrium  of  forces  are 
applied  to  thick  slices  of  the  surface  stream;  furthermore,  the  continuity 
equation  is  integrated  over  small  increments  of  time.  Simple  quadrature  for- 
mulas, notably  the  trapezoidal  rule,  are  used  to  approximate  the  integrals  by 
algebraic  expressions.  The  model  described  by  Strelkoff  and  Katopodes  (10) 
was  not  programmed  to  treat  irrigation  streams  ponded  behind  dikes  raised 
at  the  downstream  end  of  the  field.  Recession  from  the  front  of  the  irrigation 
stream,  either  before  advance  was  completed  or  after,  also  could  not  be 
modeled.  The  present  model  is  capable  of  treating  both  of  these  circumstances. 
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Additional  bases  for  nriathematical  modeling  are  provided  in  the  current 
model  by  normal-depth  kinematic-wave  theory.  If  the  depth  gradient  in  the 
surface  flow  is  sufficiently  small  compared  with  the  bottom  slope,  the  depth 
and  discharge  at  every  point  in  the  stream  are  uniquely  related  through  a 
normal-depth  stage-discharge  relationship.  This  allows  simpler,  quicker,  more 
robust^  computations  than  those  provided  by  equilibrium  (zero-inertia)  theory. 
The  range  of  applicability  of  the  kinematic-wave  model  can  be  greatly  extended 
if  application  of  kinematic-wave  theory  is  deferred  until  water  starts  to  recede 
from  the  upper  end  of  the  field.  Such  a hybrid  model  constitutes  a subset  of 
BRDRFLW,  as  does  pure  kinematic-wave  theory.  These  can  be  called  for  as 
desired  by  the  user. 

The  present  equilibrium  model  differs  also  from  its  more  recent,  undocumented 
predecessors  in  several  substantial  ways.  The  length  6x  of  cells  making  up 
the  surface  stream  has  been  divorced  from  the  time  step  6t.  For  quick,  rough 
computations  the  number  N of  cells  can  be  set  to  a small  constant,  say, 

N = 5.  Previously,  the  number  of  cells  equaled  the  number  of  time  steps  of 
advance,  which  results  in  unnecessarily  large  numbers  of  computations.  The 
present  system  also  allows  the  computational  field  length  to  be  precisely 
equal  to  the  given  field  length,  rather  than  somewhat  larger  in  accordance 
with  the  advance  achieved  in  an  integral  number  of  time  steps,  as  in  the 
earlier  models. 

Further,  irregular  bottom  profiles  can  be  input  to  the  model,  in  contrast  to  the 
strictly  plane  borders  allowed  in  the  earlier  versions.  If  an  irregular  bottom 
breaks  through  the  stream  surface,  however,  an  error  may  be  incurred,  as  the 
program  will  compute  the  motion  of  only  one  continuous  length  of  surface 
stream.  If  the  stream  breaks  up  into  several  pieces  with  islands  in  between, 
only  the  largest  piece  will  be  followed  in  detail.  The  other  pieces  are  assumed 
to  infiltrate  as  stagnant  ponds. 

The  surface-stream  profile  is  found  at  each  time  step  either  iteratively, 
through  the  Newton-Raphson  method  of  solving  the  nonlinear  governing 
equations,  or  in  a single  set  of  implicit  computations  utilizing  locally  linearized 
forms  of  the  equations.  The  latter  calculation  is  merely  the  first  step  of  the 
former.  When  used,  the  Newton-Raphson  method  typically  converges  in  about 
three  or  four  iterations. 


''Less  subject  to  computed  saw-tooth  profiles  and  spon- 
taneous aborts  of  computer  runs. 


The  skeleton  of  BRDRFLW  is  built  in  such  a way  that  variations  in  the  func- 
tional form  of  inputted  conditions  are  relatively  easily  accommodated.  For  ex- 
ample, when  in  the  course  of  computation  the  program  requires  a value  of  in- 
flow discharge  at  time  t,  a call  is  made  to  the  function  subprogram  QF(t).  Any 
desired  inflow-hydrograph  function  of  time  can  be  programmed  into  the  sub- 
program. Similarly,  the  program’s  need  for  a value  of  time  step,  say  at  time  t, 
results  in  a call  to  DTF(t,. . .);  any  desired  variation  in  time-step  size  can  be 
specified  in  that  one  subprogram.  Variable  bottom  slope  is  handled  as 
follows:  In  the  statement  of  equilibrium  among  all  forces  acting  on  the  water 
in  a cell,  the  pertinent  geometrical  feature  of  the  bottom  is  the  difference  in 
elevation  at  the  upstream  and  downstream  ends  of  the  cell.  The  subprogram 
BF(x, . . .)  provides  the  bottom  elevation  of  point  x.  Any  desired  variation  can 
be  programmed  into  that  routine.  Current  programming  allows  either  a con- 
stant bottom  slope  or  a table  of  bottom  elevations. 

The  roughness  of  the  bottom  is  characterized  by  the  Chezy  C,  which  in 
general  is  a function  of  depth  y.  The  subprogram  CHEZYC,  which  has  among 
its  arguments  y,  RDF,  a^,  and  RUFMOD,  computes  the  Chezy  C on  the  basis 
of  the  computer-supplied  y and  the  roughness  parameter  RDF,  read  in  at  the 
start  of  the  computation.  The  index  RUFMOD,  also  supplied  by  the  user, 
determines  the  nature  of  the  roughness  formula  to  be  used.  With  RUFMOD 
set  to  1,  RUF  is  interpreted  as  a constant  Chezy  C value.  With  RUFMOD  = 2, 
RUF  represents  the  Manning  n.  Some  researchers  allow  the  Manning  n to 
vary  as  a power  law  of  depth  (reflecting  the  fundamental  unsuitability  of  the 
Manning  formula  for  shallow  flow).  To  allow  comparison  with  the  results  of 
these  researchers,  the  CHEZYC  routine  utilizes  the  following  formula  when 
RUFMOD  = 2:  n = RUF  • y®";  in  the  usual  case,  a^  = 0.  With  RUFMOD  = 3, 
RUF  is  interpreted  as  the  Sayre-Albertson  x (chi)  in  their  logarithmic  formula 
(6),  considered  theoretically  more  sound  than  the  Manning  formula. 

Cumulative  infiltration  z (and  infiltration  rate)  is  also  relegated  to  a subroutine 
whose  principal  argument  is  r,  the  infiltration  time.  The  formula  currently  pro- 
grammed is  rather  general:  z = kr®  -t-  br  -i-  c,  with  a,  b,  and  c supplied  by  the 
user.  With  b = 0 and  c = 0,  the  formula  is  the  Kostiakov  function;  with 
b = 0,  the  results  represent  the  Soil  Conservation  Service  formula  (11)  typify- 
ing their  infiltration  families;  with  a = Vz,  c = 0,  the  result  is  the  Philip  for- 
mula. Other  infiltration  functions  of  time  can  be  programmed  into  this 
routine,  if  desired. 
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To  save  computation  time  in  the  ponded  case  and  with  some  irregular  bottom 
configurations,  the  computer  can  recognize  a nearly  stagnant  state  of  the  sur- 
face stream  and  compute  recession  times  directly,  based  on  a succession  of 
level  water  surfaces.  Alternately,  the  program  works  in  the  usual  mode,  com- 
puting at  successive  times  the  positions  of  the  nearly  horizontal  water  sur- 
face as  it  slowly  lowers,  as  the  result  of  infiltration,  into  the  soil. 

Input  data  can  be  entered  in  either  metric  or  English  units  or  in  dimension- 
less form.  Setting  an  index  (DMLMOD)  determines  the  system  of  nondimen- 
sionalization  used  by  the  program.  With  data  entered  in  dimensional  form,  the 
desired  set  of  characteristic  values  is  used  to  put  them  in  dimensionless 
form.  The  computations,  in  any  case,  are  carried  out  in  dimensionless  form, 
with  intermediate  results,  say,  at  each  time  step,  printed  in  dimensionless 
form.  Final  results— the  advance  and  recession  functions,  ultimate  distribu- 
tion of  infiltered  water,  runoff  volumes,  efficiencies— are  presented  in  all  sets 
of  units. 

To  enhance  the  significance  of  dimensionless  input,  the  program  computes 
the  corresponding  dimensioned  values  for  a hypothetical  border.  Hypothetical 
values  of  border  roughness,  slope,  and  inflow,  needed  to  interpret  physically 
the  given  dimensionless  input  (see  sec.  3.8),  are  entered  by  the  user,  or 
default  values  can  be  supplied  by  the  program. 

In  general,  most  input  data  defining  solution  methods  and  parameters  rather 
than  physical-problem  parameters  have  program-supplied  default  values, 
which  are  enabled  by  entering  0 (zero)  for  each  such  variable  value  called  for. 


2.  Influence  of  Computer-Facility  System  Software 


In  moving  a computer  program  from  one  facility  to  another,  certain  portions 
require  modification  because  of  system  software  peculiar  to  the  facility  that 
is  utilized  by  the  program.  BRDRFLW  was  constructed  in  FORTRAN4 
language  (MNF4  at  BKY)  using  the  computing  facilities  at  the  Lawrence 
Berkeley  Laboratory,  University  of  California  at  Berkeley  (BKY). 

The  program  was  designed  for  interactive  use,  with  a terminal  capable  of 
printing  132  characters  per  line,  such  as  the  DIGITAL  DECWRITER  III. 
Requests  for  data  are  made  in  the  form  of  prompts  to  the  operator  who  then 
enters  the  necessary  information  from  the  keyboard.  The  program  may  be  used 
in  batch  mode  by  deleting  the  subroutine  calls,  CALL  CONNECT(5)  and  CALL 
CONNECT(6)  (these  subroutines  are  part  of  the  library  at  BKY  and  allow  inter- 
active communication).  The  prompts  for  data  should  then  be  anticipated  by 
the  operator  who  submits  a small  data  deck  (typically,  a dozen  cards),  one 
card  for  each  group  of  data  requests  expected  (see  ch.  3). 

It  is  also  possible  to  allow  entry  of  input  data  from  a terminal,  while  shunting 
output  to  a line  printer.  In  the  program,  unit  number  U1  is  designated  for 
reading  input,  U2  for  printing  program  output,  and  U3  for  providing  input-entry 
prompts  to  a terminal  operator.  Numerical  values  are  assigned  to  these 
variables  in  the  data  statement  following  the  type-variable  declaration  and 
dimension  statements  in  the  main  program.  The  values  5 and  6 correspond, 
respectively,  to  input  and  output  units  at  many  computer  centers. 

In  addition,  the  program  was  designed  for  graphic  output,  in  conjunction  with 
a computer-center-library  graphics  package  (GRAFPAC  at  BKY)  and  a plotter, 
for  example,  the  TEKTRONIX  4662.  It  is  possible  to  display  depth  or  water- 
surface  profiles,  as  well  as  advance  and  recession  trajectories,  runoff  as  a 
function  of  time,  and  the  post-irrigation  infiltration  profile.  With  the  kinematic- 
wave  option  in  force,  the  wave  trajectories  in  the  x-t  plane  can  also  be  plotted. 

For  use  at  computer  centers  where  no  graphics  capabilities  exist,  reference 
to  the  GRAFPAC  routines  TVINIT,  TVWIND,  TVVIEW,  TVPLOT,  TVSEND,  and 
TVEND  must  be  deleted  from  BRDRFLW  or  dummy  routines  must  be  added, 
unless  the  software  at  the  given  center  allows  reference  in  a program  to  a 
nonexistent  routine  as  long  as  it  is  not  actually  called.  Calls  to  these 
graphics  routines  are  avoided  by  setting  all  graphics  input  parameters  (see 
Section  3.12,  LINE  12)  to  zero.  Computer  centers  with  graphics  packages  dif- 
ferent from  GRAFPAC,  will  require  calls  in  BRDRFLW  different  from  those  to 
TV. . .cited  above.  The  call  to  CCNNECT(4LFILM)  also  refers  to  BKY  software 
and  allows  interactive  graphing. 
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other  calls  in  BRDRFLW  peculiar  to  BKY  software  that  will  likely  require 
modification  for  use  at  other  computing  facilities  are  the  control  cards 
(records  1 to  9)  up  to  and  including  the  PROGRAM  card  and  DATE,  CLOCK, 
and  STATUS,  all  of  which  are  called  from  the  BRDRFLW  subroutine  RNINF. 
These  subroutines  provide  information  on  real  date  and  time  at  the  moment 
they  are  called,  the  amount  of  central  processor  time  used  up  to  that  point, 
and  other  information  pertinent  to  the  computer  run. 

To  reduce  the  difficulties  of  program  adaptation  to  various  computer  centers, 
free  format  has  been  used  but  sparingly,  specifically  only  for  input.  The  star 
(*)  in  input  statements  READ  (U1,  *). . .signifies  at  BKY  that  the  list  of  input 
data  specified  is  in  free  format,  that  is,  in  either  integer  or  decimal  form,  and 
in  fields  of  arbitrary  length,  set  off  by  commas.  This  provides  maximum  user 
flexibility  in  input  format.  All  output  format  is  achieved  through  FORMAT 
statements,  and  textual  material  therein  is  introduced  as  a specified  number 
of  Hollerith  characters. 

BRDRFLW  is  a large  program.  At  BKY  the  memory  requirement  for  compila- 
tion and  loading  is  about  176,000  words  (octal),  that  for  execution  about 
153,000  (octal). 


3.  Input 


3.1  LINE  1 (Class  of 
Computer  Run) 


As  indicated  in  chapter  2,  the  input  to  the  model  is  in  the  form  of  keyboard 
entries  during  an  interactive  session,  or  as  a series  of  punched  cards  for 
batch  processing.  The  individual  records,  call  LINES  in  the  sequel,  are  there- 
fore either  a set  of  data  preceding  a carriage-return  (send)  strike  on  a key- 
board, or  a card  of  data  input  fed  into  a card  reader.  Input  data  received  by 
the  program  are  immediately  echoed  back  to  unit  U2  (see  ch.  2 for  description 
of  input/output  units). 

Much  of  the  data  called  for  define  solution  options,  rather  than  describing 
particulars  of  a border  irrigation.  Often,  the  former  have  program-supplied 
default  values,  which  the  inexperienced  user  can  enable  by  entering  a 0 (zero) 
for  the  particular  value.  The  variables  so  endowed,  and  the  default  values  for 
each,  are  identified  in  the  sequel,  which  describes  each  variable  entry  in  turn. 
See  appendix  for  sample  input. 

The  first  line  (or  card)  of  input  is  aimed  at  interactive  users.  It  allows  the  user 
to  position  the  paper  at  the  terminal  to  whatever  vertical  placement  he  or  she 
desires  before  beginning  the  run.  For  example,  a fresh  page  can  be  brought 
up  for  recording  prompts  and  data  input,  uncluttered  by  information  pertain- 
ing to  log-in,  computer-center  messages,  and  so  on.  Then,  entry  of  a single 
digit,  between  1 and  5 inclusive,  signifies  that  the  operator  is  ready;  a 0 
brings  about  a normal  termination.  The  significance  of  this  parameter, 
TSTMOD,  follows 

TSTMOD  — indicates  the  reiationship,  if  any,  of  the  current  simulation  to 
the  last  one  executed  during  the  current  program  connec- 
tion. 

= 1 : the  current  simulation  is  independent  of  any  preceding  one; 
a complete  new  set  of  data  will  be  entered. 

= 2 : geometrical  data  are  assumed  the  same  as  in  the  preceding 
simulation.  Following  entry  of  LINE  2 (see  below),  the  pro- 
gram in  effect  jumps  to  entry  of  LINE  7.  This  option  is  con- 
venient for  management-optimization  studies.  The  manage- 
ment variables  qm,  t^o  can  be  changed  relative  to  their 
value  stored  in  the  preceding  run.  Entry  of  a zero  for  any  of 
these  causes  retention  of  the  previous  value. 

Prompts  will  be  made  for  entry  of  solution  modes  and 
parameters  (LINES  9,  10).  Zero  entries  will  cause  retention  of 
the  value  from  the  previous  run.  Diagnostic  and  plotting 
parameters  (LINES  11,  12)  must  be  entered  anew. 
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3.2  LINE  2 (Run 
Identification) 


3.3  LINE  3 (Dimension 
Control) 


= 3 : all  physical  data  are  assumed  the  same  as  in  the  preceding 
run.  Following  entry  of  the  run  identification  (LINE  2),  the 
program,  in  effect,  jumps  to  entry  of  LINE  9.  The  solution 
modes  (LINE  9)  and  parameters  (LINE  10)  can  be  changed  by 
entry  of  non-zero  values.  Zeroes  cause  retention  of  previous 
values.  Diagnostic  and  plotting  parameters  are  entered 
anew. 

= 4 : physical  parameters  and  solution  technique  are  identical  to 
those  of  the  preceding  run,  but  numerical-solution  param- 
eters (LINE  10)  can  be  changed.  Diagnostic  and  plotting 
parameters  (LINES  11,  12)  are  entered  anew. 

= 5 : all  parameters  are  retained  except  diagnostic  and  plotting 
flags  (LINES  11  and  12).  These  are  entered  anew. 

= 0 : stop 

The  program,  after  completing  an  irrigation  simulation,  or  upon  encountering 
certain  error  states,  returns  to  this  point  to  allow  entry  of  a new  set  of  condi- 
tions without  the  need  for  program  reloading.  When  the  operator  has  com- 
pleted all  the  runs  contemplated,  or  notes  his  execution-time  limit  approach- 
ing (this  information  is  supplied  by  RNINF,  prior  to  the  prompt  for  TSTMOD), 
entry  of  a 0 ends  program  control. 

The  next  line  (or  card)  of  input  requested  is  data  identifying  the  run.  This  can 
be  any  78  characters,  including  blanks,  the  user  wishes  to  have  printed  near 
the  top  of  the  printed  output  from  the  run. 

The  next  line  (or  card)  requested  contains  the  values  of  two  integer  variables, 
in  free  format,  that  is,  separated  by  commas,  defining  the  simulation  condi- 
tions. In  order,  these  are; 

INPMOD  — alerts  the  program  to  the  dimensions  of  the  irrigation  param- 
eters to  be  entered. 

= 1 ; causes  the  program  to  expect  irrigation  parameters  in  the 
metric  system:  required  depth  of  infiltration  in  cm;  field 
length  and  elevations  in  meters,  coefficient  of  power  term  in 
infiltration  formula  in  cm/hr^  in  which  a is  the  exponent, 
coefficient  of  time  term  in  infiltration  formula  in  cm/hr,  con- 
stant in  cm,  cut-off  time,  time  step,  and  maximum  irrigation 
time  computed,  in  minutes,  inflow  rate  in  liters  per  second 
per  meter  width,  flow  depths  in  meters. 
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DMLMOD 


2 : causes  the  program  to  expect  irrigation  parameters  in  the 
English  system:  required  depth  of  infiltration  in  inches,  field 
length  and  elevations  in  feet,  coefficient  of  po\wer  term  in  in- 
filtration formula  in  in/hr^,  coefficient  of  time  term  in  infiltra- 
tion formula  (final  intake  rate)  in  in/hr,  constant  in  inches, 
cut-off  time,  time  step,  and  maximum  irrigation  time  com- 
puted, in  minutes,  inflov\/  rate  in  cfs/foot  width,  flow  depths 
in  feet. 

3 : causes  the  program  to  expect  irrigation  parameters  in 
dimensionless  form. 

— The  program  works  with  dimensionless  variables.  DMLMOD 
defines  the  reference  (characteristic)  depth,  length,  time,  Yg, 
Xq,  Tg,  respectively,  computed  by  the  program  and  used, 
subsequently,  to  put  dimensioned  input  data  into  dimen- 
sionless form. 

0 : enables  one  of  the  following  two  default  values.  If  the  aver- 

age bottom  slope  of  the  border  is  greater  than  0, 
DMLMOD  is  set  to  1.  If  Sn  = 0,  DMLMOD  is  set  to  2. 

'^avg  ’ 

1 : Yo  = Yn  (normal  depth  for  the  given  inflowing  discharge  Pm, 

average  bottom  slope.  So,  and  field  roughness). 
Thus,  DMLMOD  = 1 cannot  be  used  with  horizontal 
borders. 

Xq  = yj  So 
To  = Xo/(qin/yn). 

As  a consequence,  the  dimensionless  average  bot- 
tom slope  So*  is  always  unity,  as  is  the  dimension- 
less drag  coefficient  D*.  Dimensionless  cutoff  time 
t*o  is  variable,  and  so  is  the  dimensionless  infiltra- 
tion power-law  coefficient  K*.  DMLMOD  = 1 is  the 
routine  choice  for  sloping  borders,  although 
DMLMOD  = 2 is  permissible. 

2 : Yo  satisfies  the  equation  Yq  = qfn  t^o,  in  which  Oh  is  the 

Chezy  C. 

Xq  — Qin  fco/^o 

Tq  — fco 

DMLMOD  = 2 is  the  proper  value  for  horizontal  borders. 
Then  t*o  = 1 always,  and  D*  also.  So*  is  variable  and  equals 
zero  in  horizontal  borders. 
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3.4  LINE  4 (Soil  and 
Crop  Properties) 


The  data  on  line  (card)  4 describe  the  hydraulic  properties  of  the  soil  and 
crop.  The  entries  differ,  depending  upon  whether  input  is  dimensioned  or 
dimensionless.  In  the  dimensioned  case  (INPMOD  = 1 or  2),  the  required 
variables,  eight  in  number,  are 


RUFMOD  — determines  how  program  interprets  roughness  parameters 
entered. 

= 1 : parameter  RUF  (below)  is  interpreted  as  fixed  Chezy  C. 

= 2 : parameter  RUF  is  interpreted  as  the  coefficient  in  a power 
law  of  depth  for  Manning  n;  namely,  n = RUF  • y^";  usually, 
a^  = 0 and  RUF  is  the  Manning  n. 

= 3 : parameter  RUF  is  interpreted  as  the  Sayre-Albertson  chi. 


RUF,  AN  : roughness  characteristics  of  field  (see  explanation  of 

RUFMOD,  above). 


INFMOD  — determines  how  the  program  interprets  infiltration 
parameters  entered. 

= 1 : interprets  K,  A,  B,  C (below)  as  parameters  in  the  following 
modified  power  law  for  cumulative  infiltration  z as  a function 
of  infiltration  time  r. 


z=K/+Bt+C 


A user  who  had  programmed  some  other  time  variations  of 
infiltration  would  signal  which  function  he  or  she  wished 
enabled  by  entering  an  appropriate  value  of  INFMOD. 

K,A,B,C  — infiltration  characteristics  of  field  (see  preceding  explana- 
tion of  INFMOD). 


If  dimensionless  variables  are  being  entered  (INPMOD  = 3),  the  four  infiltra- 
tion parameters  are  required,  plus  the  kind  of  roughness  formula  used;  in  the 
Sayre-Albertson  formula,  a dimensionless  x must  be  given. 

RUFMOD,  INFMOD  — same  as  above. 

K*,  A,  B*,  C*  : the  dimensionless  forms  of  K,  A,  B,  0 above. 


LINE  4a 


With  use  of  the  Sayre-Albertson  formula,  request  is  made  to  enter 

cm*  : the  dimensionless  form,  x*  = x/Yo,  of  the  Sayre-Albertson 

roughness  parameter. 
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3.5  LINE  5 (Border  The  fifth  record  describes,  in  part,  the  border  geometry.  The  three  entries  are 


Geometry) 

L 

: field  length. 

DBG 

— defines  the  downstream  boundary  condition. 

1 : open-end  border;  stream  reaching  downstream  end  of  field 
runs  off  into  a drainage  ditch. 

= 2 : blocked-end  border;  stream  reaching  field  end  is  ponded 


SOMOD 

behind  a dike. 

— defines  the  way  bottom  configuration  is  described. 

= 1 : field  is  assumed  plane;  the  bottom  slope  will  be  requested. 
= 2 : field  bottom  can  be  irregular;  pairs  of  values — distance  vs. 
elevation— will  be  requested. 

3.6  LINE(S)  6 (Bottom  Description  of  field  geometry  continues  with  specification  of  bottom  configu- 
Configuration)  ration,  in  accordance  with  SOMOD,  above.  With  SOMOD  = 1,  enter— 


SOAVG 

: the  average  slope  of  the  border. 

With  SOMOD 

= 2,  enter  instead— 

NZO 

: the  number  of  value  pairs,  distance  vs.  elevation,  that  will  be 
entered  (minimum  of  2,  maximum  of  21). 

XZ0(1),Z0(1), 
XZ0(2),  Z0(2). 
XZO(NZO), 
ZO(NZO) 

; the  NZO  value  pairs,  distance  vs.  elevation,  describing  bot- 
. .,  tom  configuration.  As  many  records  as  necessary  are  used 
to  enter  the  information.  Any  location  upstream  from  XZ0(1) 
is  assumed  to  lie  on  the  same  slope  as  between  points  1 
and  2;  any  location  downstream  from  XZO(NZO)  is  assumed 
on  the  same  slope  as  the  last  two  points. 
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3.7  LINE? 

(Management 

Parameters) 


3.8  LINE  8 (Hypothetical 
Dimensioned  Border) 


The  management  variables  are  called  for  next.  The  choice  of  variables 
depends  upon  whether  the  input  is  dimensioned  or  dimensionless.  With 
dimensional  input  (INPMOD  = 1 or  2),  the  three  entries  are — 

ZREQ  : the  required  depth  of  infiltration;  this  figure  does  not  affect 

the  advance,  recession,  runoff,  or  infiltration,  but  only  the 
assessment  of  the  merit  of  the  irrigation,  deep  percolation 
losses,  etc.  If  this  variable  is  of  no  interest^  or  it  is  desired 
to  set  the  requirement,  ex  post  facto,  to  the  minimum  or 
average-low-quarter  depths  of  infiltration,  z^eq  can  be  set  here 
simply  to  zero.  The  merit  of  the  irrigation  is  presented  in  the 
post-irrigation  synopsis  on  the  basis  of  all  three  possible 
values  of  z^gq — a given  value,  inputted  here,  the  minimum 
depth  of  infiltration,  and  the  average-low-quarter  depth  of  in- 
filtration. 

Q : volumetric  inflow  rate  per  unit  width  of  border.^ 

TCO  : time  that  inflow  is  cut  off. 

With  dimensionless  input  (INPMOD  = 3),  and  DMLMOD  = 1,  one  enters  in- 
stead, the  two  values— 

ZREQ*,  TCO*  : the  dimensionless  counterparts  of  ZREQ,  TCO. 

If  DMLMOD  = 2,  only  ZREQ*  is  called  for. 

These  records  are  requested  only  in  the  event  that  all  input  is  dimensionless 
(INPMOD  = 3).  To  increase  the  significance  to  the  user  of  entered  dimension- 
less variables,  the  variables  are  translated  by  the  program  into  corresponding 
dimensioned  variables  for  a hypothetical  border.  This  hypothetical  border  is 
specified  by  inflow,  by  roughness  if  RUFMOD  = 1 or  2,  by  slope  if 
DMLMOD  = 1,  and  by  cutoff  time,  if  DMLMOD  = 2.  These  parameters,  plus 
the  dimensionless  variables  entered  as  indicated  above  are  sufficient  to 
define  all  pertinent  dimensioned  variables  for  the  hypothetical  irrigation. 
These  are  computed  and  displayed  in  the  header  material  at  the  top  of  the 
output.  Default  values  are  supplied  for  the  required  hypothetical  values 
below,  if  zeroes  are  entered  for  the  requested  data. 


2.  . .say,  in  a nonjudgemental  study  , in  which  it  is  simply 
the  factual  result— distribution  of  infiltered  water  and 
volume  of  runoff— from  a given  combination  of  soil  and  crop 
properties,  and  design  and  management  parameters,  that  is 
needed. 
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^The  program  as  written  assumes  a constant  value  for  inflow 
until  cutoff.  Other  inflow  hydrographs  can  be  programmed 
into  the  subprogram  OF. 


INPMDH 

— determines  the  dimensional  system  in  which  hypothetical 
values  will  be  given. 

= 0 : enables  default  values  for  all  variables  in  this  group. 

= 1 : metric  system  (see  INPMOD  = 1,  sec.  3.3);  this  is  the  default 
value  of  INPMDH. 

= 2 : English  system  (see  INPMOD  = 2,  sec.  3.3). 

SOAVG 

: average  bottom  slope;  requested  only  if  DMLMOD  = 1; 
default  value  is  0.001. 

QIN 

: unit  inflow  rate;  default  value  is  6 L/sm. 

TCO 

: cutoff  time,  requested  only  if  DMLMOD  = 2;  default  value  is 
120  min. 

LINE  8a 

If  INPMDH  is  not  zero,  and  RUFMOD  is  not  3,  request  is  made  to  enter 


RUF 

: hypothetical  constant  Chezy  C,  if  RUFMOD  = 1 

or 

RUF,AN 

: hypothetical  Manning  n coefficient  and  exponent  if 

RUFMOD  = 2. 

3.9  LINE  9 (Solution-  Five  integer  solution  parameters  are  entered  next 

Mode  Parameters) 


SOLMOD 

— determines  the  type  of  solution,  see  section  3.16  for  a dis- 
cussion on  the  relative  merits  of  each  type. 

= 0 : enables  the  default  value,  SOLMOD  = 2. 

= 2 ; enables  the  zero-inertia  (equilibrium)  model. 

= 3 : enables  the  normal-depth  kinematic-wave  model. 

= 5 : enables  the  hybrid:  equilibrium  model  for  advance,  then 

when  advance  is  completed  and  recession  is  underway,  the 
kinematic-wave  model  is  used. 
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LIN  MOD  — controls  the  linearity  of  solution  at  each  time  step. 

= 0 ; enables  the  default  value,  2. 

= 1 : solution  is  locally  linearized  at  each  time  step.  This  is  the 
cheaper  and  less  accurate  option. 

= 2 ; fully  nonlinear  solution;  nonlinear  algebraic  equations  are 
solved  iteratively  at  each  time  step.  Larger  time  steps  can 
be  used.  Nonlinear  solution  is  recommended  at  least  until 
the  user  becomes  familiar  with  model  behavior.  The 
strongest  argument  for  nonlinear  solution  pertains  to  the 
zero-inertia  solution  mode.  With  LINMOD  = 2,  the  user  is 
assured  that  the  equations  of  mass  and  force  balance  are  in- 
deed satisfied  in  every  cell  at  every  time  step.  With  the 
locally  linearized  solution  (LINMOD  = 1),  one  only  knows 
that  a single  try  has  been  made  to  correct  the  first  guesses 
for  values  of  depth,  discharge,  runoff,  and  advance. 

In  the  kinematic-wave  mode,  use  of  LINMOD  = 2 ensures  treatment  of  the 
wave-trajectory  segments  as  parabolas,  rather  than  as  straight  lines. 

DTMOD  — controls  the  variation  of  the  time  increment  8t  in  the  zero- 
inertia  solution.  In  principle,  8t  can  vary.  In  the  kinematic- 
wave  solution  DTMOD  has  no  pertinence,  and  its  value  is 
ignored.  In  the  zero-inertia  solution,  whatever  the  value  of 
DTMOD,  certain  variations  in  time-step  size  occur  automati- 
cally in  the  calculations.  For  example,  if  convergence  of  a 
nonlinear,  iterative  solution  is  not  achieved  within  a specified 
number  of  iterations,  the  time  step  is  automatically  cut  in 
half.  (Also,  see  ch.  6,  pt.  II.) 

= 0 : enables  the  default  value,  2. 

= 1 : 6t  = constant  = DT(STD),  entered  in  line  10. 

= 2 : DTfSTD)  is  divided  into  N(STD)''  equal  parts  to  form  a smaller 
6t.  The  small  increment  is  used  for  the  first  N(STD)  time 
steps.  During  this  period,  the  number  of  cells  comprising  the 
stream  increases  at  each  step,  from  one  to  N(STD).  Sub- 
sequently, the  number  of  cells  remains  constant,  while  the 
time  step  gradually  increases  until  the  full  value  DT(STD)  is 
achieved.  Thenceforth,  except  as  noted  above,  this  serves  as 
the  time  step.  This  procedure  is  designed  to  minimize  the  in- 
fluence of  the  assumed  shape  of  the  stream  profile  in  the 
first  time  step. 


. . .entered  in  line  10. 


ISUPZA  — Limitations  in  the  numerical  solution  can  lead  to  the  exclu- 
sion of  portions  of  the  surface  stream  from  computation.^  In 
the  event  islands  are  formed,  BRDRFLW  chooses  the  longest 
length  of  stream  on  which  to  continue  the  calculations.  The 
remainder  of  the  stream  is  excluded.  The  excluded  portions 
can  be  assumed^  to  infiltrate  approximately  as  stagnant 
ponds,  to  preserve  an  overall  mass  balance. 

= 0 : the  default  value,  causes  such  an  adjustment  of  the  reces- 
sion curve  as  to  effect  at  each  excluded  node  the  addition 
of  the  surface  depth  there  to  the  local  infiltrated  depth.  The 
user  is  cautioned  that  this  procedure  may  lead  to  irregular 
aberrations  in  the  recession  curve,  say  if  the  excluded  por- 
tions of  the  stream  profile  exhibit  large-amplitude  saw  teeth. 

= 1 : suppresses  these  additions  to  the  infiltration  profile.  The 
consequent  neglect  of  some  of  the  stream  volume  leads  to 
volume  error,  as  will  be  noted  in  the  program  output. 

ZADMOD  — controls  the  nature  of  the  solution  after  stream  flow  has 

become  essentially  stagnant.  With  very  low  infiltration  rates, 
and  little  or  no  runoff,  very  great  irrigation  times  can  result. 
Array  dimensions  in  the  program  limit  the  number  of  time 
steps  to  800. 

= 0 : enables  the  default  value,  2. 

= 1 : normal  procedure  is  followed  regardless  of  the  values  of 
flow  velocities  encountered. 

= 2 : the  program  monitors  the  water  surface  in  the  stream.  When 
this  becomes  essentially  level,  recession  is  assumed  to 
occur  as  for  a still  pond  (or  a series  of  ponds,  if  the  bottom 
is  sufficiently  irregular).  The  water  surface  in  each  pond  is 
assumed  to  fall  at  a time-varying  rate  equal  to  the  current 
distance-averaged  infiltration  rate  in  that  pond. 


^Islands  can  form  in  shallow  flow  over  a border  with  an 
irregular  (nonplane)  bottom.  Further,  even  in  a plane  border, 
an  imperfect  calculation  can  lead  to  saw-tooth  stream  pro- 
files. Computation  of  a negative  depth  anywhere  implies  for- 
mation of  an  island.  BRDRFLW  is  designed  to  compute  only 
a single  length  of  flowing  stream;  surges  in  tandem, 
separated  by  islands,  are  not  allowed. 

. .with  varying  degree  of  error. 
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3.10  LINE  10  (Numerical 
Solution  Parameters) 


Numerical  parameters  governing  the  behavior  of  the  solution  are  entered 
next,  four  with  the  locally  linearized  solution,  five  with  the  fully  nonlinear 
solution. 

N(STD)  : the  number  of  increments  (slices)  into  which  the  stream 

length  is  divided.  Array  dimensions  limit  this  value  to  80. 
Typical  values  are  10  to  20.  A value  of  5 gives  relatively 
crude,  inexpensive  results.  0 (zero)  enables  the  default  value 
of  20. 

RDX  : the  ratio  of  adjacent-cell  lengths.  A value  less  than  unity 

leads  to  a gradual  decrease  in  cell  length,  from  the  middle 
of  the  stream  outwards  towards  the  stream  boundaries. 
Small  cell  lengths  in  regions  of  sharp  profile  curvature  make 
the  choice  of  </>  weighting  factors  (see  sec.  6.2)  less  critical 
than  it  would  be  otherwise.  Thus,  small  cells  can  be  concen- 
trated in  regions  where  they  are  most  needed. 

The  value,  unity,  leads  to  uniform  cell  size.  A zero  value 
entered  for  RDX  leads  to  the  default  value  0.8,  if  N(STD)  is 
less  than  20,  and  0.9  if  N(STD)  is  greater  than  20. 

Too  small  a value  of  RDX  can  lead  to  a cell  length  so  small 
at  the  stream  front,  that  during  advance  it  is  exceeded  by 
the  correction  to  the  first-guess  advance  increment.  When 
this  happens,  the  nodes  are  automatically  repositioned 
within  the  corrected  stream  length  and  the  calculation  of  the 
time  step  is  begun  anew. 

Furthermore,  a small  value  of  NSTD  coupled  with  small  RDX 
leads  to  large  cell  lengths  in  the  interior  of  the  flow.  This 
can  lead  to  inaccuracies  if  the  profile  is  not  smooth  and 
gradually  changing  in  the  interior.  For  example,  in  the  event 
of  ponding  behind  a dike  blocking  the  end  of  a relatively 
steep  border,  the  juncture  between  the  fast  moving  stream 
in  the  upstream  portions  and  the  nearly  stagnant  pool  in  the 
downstream  portion  can  make  a relatively  sharp  angle, 
rather  than  a smooth  and  gradual  transition. 

DT(STD)  : standard  time  step  upon  which  5t  is  based  in  accordance 

with  DTMOD,  entered  in  line  9. 


16 


TMAX 


J.11  LINE  11 
Diagnostics) 


: The  computation  will  stop  when  computed  irrigation  time  ex- 
ceeds TMAX.  Such  an  upper  limit  is  introduced  to  allow  a 
limited  study.  TMAX  also  prevents  excessive  computation  in 
case  a peculiar  combination  of  input  parameters  is  entered 
(for  example,  ponded  border  with  nearly  zero  infiltration 
rate).  To  avoid  exceeding  array  dimensions,  which  allow  a 
maximum  of  800  time  steps,  set  TMAX  < 800  • DT(STD) 
(when  time  step  is  constant).  A 0 (zero)  enables  the  default 
value,  approximately  800  • DT(STD). 

If  the  fully  nonlinear,  iterative  solution  (LINMOD  = 2)  is  used,  the  number  of 
iterations  must  be  limited.  The  user  supplies  this  maximum  number,  namely— 

UMAX  : (requested  only  if  LINMOD  = 2);  limits  the  number  of  itera- 

tions in  selected  nonlinear  procedures;  rarely  are  more  than 
10  required.  20  is  a recommended  value  to  allow  the  pro- 
gram more  leeway  in  convergence  problems.  0 (zero)  enables 
the  default  value,  JMAX  = 20. 

Four  integer  entries  and  one  real  entry  are  required  in  this  record,  which  con- 
trols the  output  of  diagnostic  information. 

IDIAG  — determines  the  level  of  diagnostic  information  printed,  to 

facilitate  troubleshooting.  A given  level  enables  all  lower 
levels  as  well. 

= 0 : almost  no  diagnostic  information  is  printed  — merely  the  in- 
put in  dimensional  and  dimensionless  form,  advance  and 
recession  data,  and  ultimate  distribution  of  infiltered  water, 
efficiencies,  and  so  forth.  Evidences  of  irregularities  in  the 
solution,  however,  are  also  printed,  to  alert  the  user. 

= 1 : a brief  survey  of  the  stream  at  each  time  step  (or  for  each 
kinematic-wave  trajectory)  is  printed.  When  the  computation 
of  the  irrigation  has  concluded,  also  printed  is  a listing  of 
the  advance  and  recession  trajectories  as  computed  at  each 
time  step. 

= 2:  Zero-inertia  mode 

. . .informs  the  user  that  the  time  step  has  been  automati- 
cally cut  in  half,  following  computational  difficulties. 

Kinematic-wave  mode 

. . .enables  print  of  runoff  calculations  for  each  kinematic- 
wave. 

. . .in  the  hybrid  mode,  enables  print  of  profile  used  as  start- 
ing point  for  kinematic-wave  analysis;  also  prints  advance 
curve  as  computed  by  zero-inertia  model. 
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= 3 : Zero-inertia  mode 

. . .enables  printing  of  a stream  profile  at  each  time  step. 

Kinematic-wave  mode 

. . .enables  printing  of  variation  of  stream  variables  along 
the  kinematic-wave  trajectory,  for  each  trajectory  calculated. 

= 4;  Zero-inertia  mode 

. . .in  the  ponded  case,  shows  computations  leading  to  first 
guess  for  downstream  pool  elevation  in  the  first  step  after 
stream  has  reached  field  end. 

Kinematic-wave  mode 

. . .in  the  ponded  case,  shows  computations  leading  to  pool 
dimensions  at  start  of  standing-pool  recession. 

Higher  levels  of  diagnostics  allow  the  user  to  explore  the  internal  workings  of 
the  program  in  ever  greater  detail.  Level  5 refers  both  to  the  zero-inertia  and 
kinematic-wave  modes  of  solution,  while  values  of  IDIAG  equal  to  6,  7,  or  8 
refer  only  to  the  zero-inertia  mode. 

Kinematic-wave  mode 

= 5 : prints  the  results  of  each  iteration. 

Zero-inertia  mode 

= 5 : prints  the  results  of  each  iteration  in  selected  computations. 
In  particular,  at  each  iteration,  it  prints  residuals  from  both 
continuity  and  force-equilibrium  equations,  as  well  as  values 
of  stream  variables  and  predicted  and  resulting  values  of 
pressure,  weight,  and  drag  in  the  force  equation.  Because  of 
its  high  degree  of  nonlinearity,  the  force  equation  usually  ex- 
hibits the  greatest  lack  of  correspondence  between  pre- 
dicted and  actual  variable  values.  This  level  of  diagnostic 
produces  one  line  of  data  for  each  cell. 

Level  5 also  enables  printing  of  the  new  cell  dimensions, 
after  the  computational  stream  length  has  been  shortened. 

= 6 : prints  the  six  components  of  the  continuity  and  force  equa- 
tions, the  sensitivity  of  each  to  each  of  the  stream  variables, 
and  the  incremental  changes  to  these  variables,  for  each 
cell,  in  each  iteration  (six  lines  of  data  for  each  cell). 

= 7 ; prints  the  raw  output  from  the  double-sweep  routine. 

= 8 : displays  the  coefficients,  solution  vector,  and  right-hand 
sides  of  the  set  of  linear  equations  solved  at  each  iteration. 


IDCH 

: the  number  1 of  the  time  step  at  which  it  is  desired  to 

change  IDIAG;  for  example,  a program  may  work  well  up  to 
the  132nd  time  step,  then  some  irregularities  are  noted. 

IDIAG  may  be  held  to,  say,  zero  or  one  for  the  first  132  time 
steps,  and  IDCH  set  to  132  to  cause  a change  in  the  value  of 
IDIAG  at  the  132nd  step;  if  no  changes  in  IDIAG  are  desired, 
set  IDCH  = 0. 

ID2 

; the  new  value  of  IDIAG,  from  1 = IDCH  onward.  If  no 
changes  in  IDIAG  are  desired,  set  ID2  = 0. 

IPRZA 

— for  certain  evaluations  of  the  ultimate  profile  of  infiltrated 
depths,  it  Is  desirable  to  construct  an  equivalent  profile  with 
monotonic  increasing  depths.  The  program  performs  this 
reordering. 

= 1 : enables  printing  of  equivalent  profile. 

= 0 : suppresses  such  printing. 

FLGPVE 

— value  of  overall  relative  volume  error,  which  if  exceeded  will 
cause  a special  print  message. 

= 0 : the  default  signal,  which  sets  FLGPVE  = 0.015. 

.12  LINE  12  (Plotting  This  line  (card)  of  data  contains  5 integers,  all  related  to  plotting.  If  no  plots 
'arameters)  are  desired,  or  in  the  absence  of  facilities  for  automatic  plotting,  all  5 items 

should  be  set  to  zero.  Descriptions  of  these  variables  follow: 


IPLOTW  : 

= 1 : advance,  recession,  and  runoff  data  are  plotted  as  functions 
of  time;  ultimate,  post-irrigation  infiltration  depths  are  plotted 
vs.  distance;  = 0:  no  plotting  of  these  variables. 

IPLOTY  : 

= N : with  positive  integer  N,  depth  or  water-surface-elevation  pro- 
files are  plotted  at  intervals  N;  =0:  no  plots  of  depth  or 
water  surface. 

IPLOTH  : 

= 1 : water-surface  elevation  rather  than  depth  is  plotted;  = 0: 
water  surface  not  plotted. 

IPLOTC  : 

= 1 : trajectories  of  computed  kinematic  waves  are  plotted;  = 0: 
no  plots  of  these  waves. 
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IPWAIT 

— sometimes,  in  interactive  use,  a pause  is  needed  before 
automatic  plotting  of  each  profile  begins;  this  allows  some 
adjustment  in  the  plotter,  for  example,  a change  in  pen 
color. 

= 1 : prompts  the  user  to  enter  a zero,  before  plotting  com- 
mences. Following  any  necessary  adjustments,  the  user 
enters  0,  and  the  profile  is  plotted.  Prior  to  plotting  the  next 
profile,  the  prompt  appears  again. 

= 0 : profiles  are  plotted  one  after  the  other,  at  intervals  dictated 
only  by  the  rate  of  data  transmission. 

3.13  LINE  13a  (Plot 

Scales) 

FSX 

(Called  for  only  if  profile  plotting  is  enabled) 

: full-scale  value  of  abscissa;  a zero  enables  a default  value. 

FSY 

: full-scale  value  of  ordinate;  a zero  enables  a default  value. 

LINE  13b 


(called  for  only  if  trajectory  plotting  is  enabled) 

PLTXMX 

: full-scale  value  of  the  x-coordinate.  A zero  here  enables  the 
default  value,  border  length.  A negative  entry  allows  reten- 
tion of  value  from  preceding  run. 

PLTTMX 

: full-scale  value  of  the  time  coordinate.  A zero  enables  the 
default  value.  A negative  entry  allows  retention  of  value  from 
preceding  run. 

3.14  LIN E(S)  14  (Dummy  MP 

Pause  Variable) 

: enter  0 repeatedly,  before  every  profile  plot  if  IPWAIT  = 1 
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.15  Summary  of  input 


Record  Parameter 


Function 


Permissible 

values 


Defaults 


Line  1: 

Class  of  computer  run 

TSTMOD 

terminate,  or  continue  with  0,1, 2,3,4, 5 

certain  data  stored,  some  new 
data  entered 

Line  2: 

Run  identification  any  78  blanks 

characters 

Line  3: 

Dimension  control 

INPMOD 

DMLMOD 

indicates  units  of  input  data  1,2,3  — 

controls  form  of  dimensionless 

variables  0,1,2  0 — 1 or  2 

Line  4: 

Soil  and  crop  properties 

With  dimensioned  input  (INPMOD  = 1,2): 

RUFMOD  determines  roughness  formula  1,2,3 


RUF 

AN 

INFMOD 

K 

A 

B 

0 

roughness  coefficient  user  entry 

exponent  for  variable  Manning  n user  entry  0 

determines  infiltration  formula  1 

I terms  in  z = Kr^  + Br  + C user  entry 

or:  With  dimensionless  input  (INPMOD  = 3): 

RUFMOD  determines  roughness  formula  1,2,3 


INFMOD 
K*  ) 

B*  [ 

C*  ' 

determines  infiltration  formula  1 

terms  in  z*  = K*t*^  + B*t*  + C* 

Line  4a:  (called  for  only  if  INPMOD  = 3,  RUFMOD  = 3) 

CHI*  dimensionless  Sayre-Albertson  X user  entry 


Line  5: 

Border  geometry 

L 

DBG 

SOMOD 

field  length  user  entry 

downstream  boundary  condition  1,2 

plane,  or  irregular  bottom  1,2 
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Record  Parameter 


Function 


Permissible 

values 


Defaults 


Line  6:  Bottom  configuration 

For  plane  bottom  (SOMOD  = 1): 

With  dimensioned  input  (INPMOD  = 1,2): 

SOAVG  bottom  siope  user  entry 


or:  With  dimensionless  input  (INPMOD  = 3,  DMLMOD  = 2): 

SO*  dimensionless  bottom  slope  user  entry 


or:  For  irregular  bottom  (SOMOD  = 2): 

NZO  number  of  value  pairs  (x^^,  Zq) 
to  be  entered 

XZ0(1),  Z0(1),  XZ0(2),  Z0(2),.  . XZO(NZO),  ZO(NZO) 


Line  7: 


Management  parameters 


With  dimensioned  input  (INPMOD  = 1,2): 

ZREQ  required  depth  of  infiltration 
Q unit  inflow  rate 

TOO  cutoff  time 

or:  with  dimensionless  input  (INPMOD  = 3) 

. . .and  DMLMOD  = 1: 

ZREO*  dimensionless  required  depth 
TOO*  dimensionless  cutoff  time 


2-21 

(inclusive) 
user  entry 


user  entry 
user  entry 
user  entry 


user  entry 
user  entry 


or:  . . .and  DMLMOD  = 2 

ZREO*  dimensionless  required  depth  user  entry 

Line  8:  Hypothetical  dimensioned  border  (called  for  only  if  INPMOD 


= 3) 


With  DMLMOD  = 1 


INPMDH 

controls  units  of 
hypothetical  example 

0,1,2 

0 - 1 

SOAVG 

hypothetical  average 
bottom  slope 

> 0 

0 - 0.001 

OIN 

hypothetical  unit  inflow  rate 

0-6  L/sm 

or:  With  DMLMOD 

INPMDH 

= 2 

controls  units  of 
hypothetical  values 

0,1,2 

0 - 1 

Permissible 


Record 

Parameter 

Function 

values 

Defaults 

QIN 

hypothetical  unit  inflow  rate 

user  entry 

0—6  L/sm 

TOO 

hypothetical  cutoff  time 

user  entry 

0 — 120  min 

Line  8a: 

(called  for  only  if  INPMDH  ^ 0,  RUFMOD  ^ 3) 

With  RUFMOD 

= 1 

RUF 

hypothetical  constant  Chezy  0 

user  entry 

0 - 12.55 

or:  With  RUFMOD 

= 2 

RUF 

hypothetical  Manning-n 

coefficient 

user  entry 

0 - 0.05 

AN 

hypothetical  Manning-n 

exponent 

user  entry 

0-0 

Line  9: 

Solution-mode  parameters 

SOLMOD 

solution  technique 

0,2, 3, 5 

0-2 

LINMOD 

solution  linearity 

0,1,2 

0-2 

DTMOD 

controls  time-step  variation 

0,1,2 

0-2 

ISUPZA 

controls  option  for  treatment  of 

excluded  portions  of  stream 

0,1 

0 

ZADMOD 

controls  option  for  stagnant 

ponding 

0,1,2 

0-2 

Line  10: 

Numerical  solution  parameters 

N(STD) 

number  of  stream  increments 

0-80 

0-20 

RDX 

spatial  rate  of  cell-size 

00 

o 

1 

o 

decrease 

0-1 

or  0.9 

DT(STD) 

basic  time  step 

user  entry 

TMAX 

maximum  irrigation  time 

0 - 800  • 

allowed 

user  entry 

DT(STD) 

UMAX 

(requested  only  if  LINMOD  = 2) 

maximum  number  of  iterations 

user  entry 

o 

CM 

! 

o 

Line  11: 

Diagnostics 

IDIAG 

controls  amount  of  diagnostic 

0 - 

information 

0-8 

summary 

only 

IDCH 

time  step  at  which  IDIAG  is 

changed 

user  entry 

0 - 

no  change 
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Record  Parameter 


Function 


Permissible 

values 


Defaults 


IDIAG2 

new  level  of  diagnostics  after 

change 

0-8 

0 

IPRZA 

if  ^ 0,  prints  equivalent 
monotonic  increasing  ultimate- 

infiltrated-depth  distribution 

0,1 

0 — no  print 

FLGPVE 

flags  overall  relative  volume 
error  in  excess  of  FLGPVE 

> 0 

0 

1 

o 

o 

L 

cn 

Plotting  parameters  (total  of  5) 

IPLOTW 

controls  plot  of  advance, 
recession  trajectories,  runoff 
function,  ultimate  infiltrated- 

0,1 

0:no  plot 

depth  profile 

1:plot 

IPLOTY 

controls  frequency  of 

integer 

0:no  plot  of 

plotting  depth  or  water-surface 

depth 

elevation 

profile 

IPLOTH 

controls  plot  of  water-surface 

0,1 

0:no  plot 

elevation  profile 

1:plot 

IPLOTC 

controls  plotting  of  kinematic- 

wave  trajectories 

0,1 

// 

IPWAIT 

allows  a pause  before 

0,1 

0:no  pause 

profile  plot 

lipause  to 

read  a 
dummy 
variable 

Profile-plot  scales 

FSX 

full-scale  value  of  abscissa 

user  entry 

0 

FSY 

full-scale  ordinate  value 

user  entry 

0 

Line  13b:  Trajectory  plot  scales 

PLTXMX  full-scale  x-coordinate  value  -1,or 

user  entry 

PLTTMX  full-scale  t-coordinate  value  - 1,  or 

user  entry 

( - 1 leads  to  retention  of  values  stored  in  a previous  run) 

Line  14:  Dummy  pause  variable  (called  for  only  if  IPWAIT  = 1) 

MP  entered  repeatedly,  before 

every  profile  plot 
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3.16  Relative  Merits  of 
the  Submodels  (With 
Reference  to  SOLMOD, 
Section  3.9) 


For  the  interactive  user,  the  program  contains  a bail-out  feature  that  can  be 
used  if  incorrect  data  have  been  inadvertently  transmitted,  and  the  mistake  is 
noted  before  completion  of  data  input.  Entry  of  - 1 for  the  first  variable  value 
(and  appropriate  numbers  for  the  remaining  entries  in  the  line)  in  any  of  lines 
3,  4,  4a,  5,  7,  8,  8a,  9,  10,  11,  12  returns  the  user  to  the  beginning  of  data 
input,  without  loss  of  program  control. 

The  zero-inertia  technique  is  the  usual  for  computations  of  border-irrigation 
streamflow.  Under  conditions  of  very  steep  bottom  slope  (typically  greater 
than  about  one-half  pet.),  poor  numerical  approximations  in  the  solution 
algorithms  (for  example,  determination  of  average  value  of  drag  on  a cell  of 
shallow  water  from  drag-per-unit-length  values  at  the  upstream  and  down- 
stream boundaries  of  the  cell,  or  determining  the  weight  of  the  most  down- 
stream cell  from  the  depth  at  its  upstream  end,  its  length,  and  an  assumed 
shape  factor  0y)  lead  to  computed  profiles  exhibiting  a saw-tooth  character.  If 
the  amplitude  of  the  oscillations  is  not  too  great  a fraction  of  the  depth,  this 
does  not  constitute  a serious  problem,  because  the  computation  is  funda- 
mentally stable,  and  the  irregularities  do  not  tend  to  grow  at  each  time  step. 
However,  if  a computed  depth  drops  anomalously  below  zero,  somewhere  in 
the  interior  of  the  stream,  the  computed  recession  will  occur  too  quickly,  and 
the  relative  volume  error  of  the  computation  will  become  high  (say,  substan- 
tially greater  than  -i-O.OI).  While  such  problems  can  often  be  alleviated  by 
the  use  of  smaller  time  steps  and  cell  lengths,  the  cost  of  computations  then 
increases.  An  alternate  approach  to  the  problem,  since  it  arises  with  high 
values  of  bottom  slope,  utilizes  the  concepts  of  kinematic-wave  theory. 

Kinematic-wave  theory  based  on  a normal-depth  stage-discharge  relationship 
is  an  inexpensive,  robust  technique  and  yields  results  very  close  to  those  of 
zero-inertia  theory,  provided  the  depth  gradient  is  everywhere  small  compared 
with  the  bottom  slope.  This  eliminates  kinematic-wave  theory  completely  in 
horizontal  borders.  But  in  sloping  borders,  during  advance,  the  upstream 
depth  is  of  the  order  of  normal  depth,  so  that  by  the  time  the  stream  has 
advanced  to  a distance,  say,  20  times  yJSo,  the  depth  gradient  is  less  than 
5 percent  of  the  bottom  slope.  In  practice,  once  L*  = L So/yn  is  greater  than 
about  10,  advance  is  accurately  computed  by  kinematic  wave  theory. 

Right  after  inflow  cutoff,  however,  the  upstream  depth  drops  sharply  while 
depths  in  the  interior  of  the  stream  reduce  more  gradually.  Thus,  there  is  a 
substantial  depth  gradient  in  the  stream  for  a while  after  cutoff.  This  leads  to 
inaccurate  computation  of  recession  by  kinematic-wave  theory,  even  at 
substantial  values  of  L*. 
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If,  however,  application  of  kinematic-wave  theory  is  delayed  until  recession  as 
computed  by  zero-inertia  theory  begins,  the  depth  gradient  has  generally  had 
the  opportunity  to  decrease  in  magnitude  to  the  point  where  kinematic-wave 
theory  becomes  applicable.  This  is  the  basis  of  the  hybrid  model.  In  operation 
of  this  model,  application  of  kinematic-wave  theory  is  delayed  until  both 
recession  starts  and  the  advancing  stream  has  reached  field  end.  This  simpli- 
fies the  computation  of  the  kinematic-wave  trajectories,  which  otherwise 
would  have  to  be  computed  simultaneously  with  the  trajectory  of  the  advanc- 
ing front.  The  hybrid  model  exhibits  good  accuracy  for  L*  as  low  as  unity. 

As  a comparison  of  typical  execution  costs,  if  the  zero-inertia  model  costs 
$1.00  to  run,  the  kinematic-wave  model  could  be  run  for  $0.16  and  the  hybrid 
for  $0.36. 


4.  Printed  Output 


The  first  output  to  appear  following  standard  printouts  originating  with  the 
computer  system  will  be  the  line.  . ***BORDER-IRRIGATION-FLOW 
PROGRAM***.  . This  is  followed  (provided  the  subprogram  RNINF  or  its 
counterpart  in  the  given  computer  system— see  ch.  2— is  in  operation)  by  cer- 
tain real,  on-line  computing  information:  calendar  date  and  time,  elapsed  cen- 
tral processor  time,  and  computation  units  used  and  remaining  available  for 
use.  This  is  followed  by  prompts  for  the  input  data  and  the  data  themselves, 
if  these  have  been  routed  to  the  output  printing  unit  (ch.  See  appendix  for 
sample  printout. 

Formal  output  begins  with  the  line: 

*********************  border  irrigation  flow  ********************* 

Then,  if  the  subprogram  RNINF  or  its  counterpart  is  in  operation,  the  calendar 
date  and  time  are  printed. 

This  is  followed  by  the  run  identification  entered  as  data  in  LINE  2. 

The  header  information  that  follows  is  divided  into  five  groups.  The  first  three 
describe  the  problem;  the  fourth  lists  the  computed  characteristic  reference 
parameters  which  enter  into  the  dimensionless  problem  parameters,  and  the 
last  group  describes  the  parameters  of  the  numerical  solution. 

The  problem  parameters  consist  of  (1)  the  hydraulic  properties  of  crop  and 
soil,  that  is,  roughness  and  infiltration,  (2)  the  border  geometry,  that  is,  bot- 
tom configuration,  and  (3)  the  management  parameters,  that  is,  required 
depth  of  application,  inflow  rate,  and  cutoff  time. 

The  characteristic,  reference-parameter  information  consists,  first,  of  a state- 
ment (governed  by  DMLMOD— see  sec.  3.3)  defining  the  reference  variables. 
Next,  normal  depth  for  the  inflow,  average  slope,  and  roughness  and  Froude 
number  at  normal  depth  are  printed  (if  the  slope  is  greater  than  zero).  These 
are  followed  by  the  characteristic  depth,  length,  discharge,  and  time. 

The  parameters  XOP,  TOP,  P refer  to  the  system  of  nondimensionalization 
used  in  (5),  and  are  not  used  in  the  program;  they  are  included  only  for  com- 
parison with  the  earlier  work. 

The  solution  parameters  define  the  solution  technique,  and  also  the  numeri- 
cal elements  thereof;  namely,  the  number  of  cells,  size  of  time  step,  and  max- 
imum number  of  iterations  allowed  in  a nonlinear  solution. 


^In  the  event  DMLMOD  = 2 and  IDIAG  > 5,  the  progress  of 
the  iterations  for  the  characteristic  depth  Yq  will  be  printed 
next;  this  will  not  usually  be  the  case. 
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In  general,  the  first  column  of  data  appears  in  English  units,  the  second  is  in 
metric  units;  dimensioniess  forms  appear  in  the  third  coiumn.  in  the  event 
that  data  are  entered  in  dimensioniess  form,  corresponding  dimensioned 
values  computed  for  a hypothetical  border  irrigation  are  also  printed  for 
added  significance. 


The  next  series  of  output  records  to  appear  depends 
for  IDIAG,  IDCH,  and  ID2  (see  sec.  3.11). 


upon  the  values  entered 


In  the  output  for  IDIAG  = 1,  the  index  J refers  to  the  number  of  iterations,  I is 
the  number  of  time  steps  (t  = 0 when  i = l),  LB  refers  to  the  ieft  (upstream) 
boundary  of  the  stream,  and  RB  refers  to  the  right  (downstream)  boundary  of 
the  stream.  FIYLN  and  FIYRN  refer  to  the  shape  factors  0y  appiied  at  the  ieft 
and  right  boundaries,  respectiveiy,  of  the  iast  computationai  ceii.  The  variabie 
TOOTH  is  a measure  of  the  irregularity,  or  iack  of  smoothness,  of  the  com- 
puted surface-water  profile.  It  represents  the  largest  value,  in  the  absolute 
sense,  of  the  variable 


s - (yk  + 2 + yk  - 2yk^i)  - (y,^^  + _ 2yj  (4 -I) 

within  the  profile.  KTOOTH  is  the  value  of  k for  which  s is  maximum.  The 
variable  PCVERR  is  the  relative  volume  error  incurred  over  the  extant  period 
of  irrigation;  it  is  the  volume  error  divided  by  the  volume  of  inflow  The 
volume  error  is  the  inflow  volume,  minus  the  sum  of  surface,  infiltrated  and 
runoff  volumes.  All  diagnostic  data  are  given  in  dimensionless  form.  If  needed 
dimensional  values  can  be  obtained  through  multiplication  by  the  appropriate 
characteristic  parameter  (QO,  YO,  XO,  or  TO)  listed  in  the  header  section  of 


At  the  conclusion  of  the  irrigation-modeling  computations,  the  results  of  post- 
irrigation  analysis  (see  sec.  5.5)  are  printed— the  ultimate  infiltrated-depth 

profile,  measures  of  the  merit  of  the  computation,  and  a synopsis  of  the  irri- 
gation results. 

As  a general  note— appearance  of  the  code  number  -7777  in  place  of  a 
variable  value  signifies  that  the  output  variable  in  question  has  not  been 

properly  defined,  for  example,  if  the  definition  formula  contains  a divide  bv 
zero. 


The  final  bit  of  output  from  the  run  is  a statement  of  the  number  of  comput- 
ing units  used  for  the  run  (provided  subprogram  RNINF  or  a substitute  is 
attached).  The  program  then  returns  to  the  beginning  with  a prompt  for 
TSTMOD  (see  sec.  3.1). 
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5.  Overview  of  Instruction  Sequence  in  BRDRFLW: 
Major  Branches 


5.1  Input 


5.2  SOLMOD  = 2 
Zero-Inertia 
(Equilibrium)  Model 


5.3  SOLMOD  = 3 

Kinematic-Wave 

Analysis 


Input  data  is  requested  and  read.  Data  entered  in  the  English  system  is  con- 
verted to  metric  (SI).  Characteristic  variables  are  defined,  and  input  data  are 
converted,  if  necessary,  to  dimensionless  form.  Dimensioned  input  variables, 
computed  characteristic  variables,  and  dimensionless  input  variables  are 
printed. 

The  program  then  branches  on  the  choice  of  solution  mode,  SOLMOD. 

For  each  of  a succession  of  time  steps,  the  profile  of  depth  and  discharge  in 
the  surface  stream  is  computed.  Each  profile  except  the  first  is  based  on  the 
results  of  the  preceding  time  step.  The  computations  are  based  on  satisfac- 
tion of  mass  conservation  and  equilibrium  of  forces  in  each  of  the  N slices 
(cells)  comprising,  together,  the  total  length  of  the  surface  stream.  The  pro- 
files gradually  lengthen  until  the  stream  front  reaches  the  end  of  the  field 
(under  certain  input  conditions,  advance  will  halt  spontaneously,  short  of  the 
field  end).  Stream  length,  more  specifically  the  location  of  the  stream  front  at 
each  instant  of  time,  is  stored  for  subsequent  construction  of  the  advance 
curve.  Once  the  end  of  the  field  is  reached,  either  runoff  or  ponding  takes 
place  there;  runoff  rate  is  monitored  and  integrated  over  time  to  yield  the 
volume  of  runoff.  Some  time  after  cut  off  of  the  inflowing  stream,  the  com- 
puted profiles  begin  to  shorten,  with  recession  occurring  from  either  the 
upstream  or  downstream  end  of  the  stream,  or  both.  Locations  of  the  stream 
boundaries  are  stored  for  construction,  ultimately,  of  the  recession  curve. 

With  IDIAG  = 1,  a line  of  output  descriptive  of  the  stream  is  printed  at  every 
computational  time  step.  With  IDIAG  = 3,  the  computed  profiles  of  dimen- 
sionless bottom  elevation  b,  depth  y,  water  surface  elevation  h,  discharge  q, 
and  depth  of  infiltration  z are  also  printed.  When  the  stream  length  has  been 
reduced  to  zero,  the  irrigation  is  over,  and  the  program  skips  to  Post-Irrigation 
Analysis  (sec.  5.5).  The  zero-inertia  computations  are  described  in  detail  in 
chapter  6. 

As  described  in  chapter  7 (also  Bassett  and  others  (1)  and  the  pioneering 
paper  (7)),  kinematic-wave  theory  applied  to  the  partial  differential  equation  of 
mass  conservation  (equation  of  continuity)  and  a given  stage-discharge  rela- 
tion (say,  normal-depth  assumption  everywhere  in  the  stream)  yields  a simul- 
taneous pair  of  ordinary  differential  equations.  One  describes  a curve  in  the 
x-t  plane— the  so-called  kinematic-wave  trajectory— in  terms  of  the  depth  at 
points  on  the  curve,  while  the  other  describes  the  variation  in  depth  along  the 
curve.  During  advance,  the  kinematic-wave  trajectories  end  abruptly  when 
they  reach  the  trajectory  of  the  wave  front,  with  an  abrupt,  step  rise  in  depth 
and  discharge  from  zero  values  on  the  field  ahead  of  the  wave  to  non-zero 
values  at  the  wave  front,  in  the  wave.  At  the  wave  front,  water  velocity  in  the 
wave  and  speed  of  propagation  of  the  wave  front  down  the  border  are  equal. 
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Kinematic-wave  trajectories  are  constructed  in  the  x-t  plane  separated  by  con- 
stant time  intervals  at  the  left  boundary  x = 0,  from  whence  they  emanate; 
depths  (and  discharge)  are  known  there.  Their  intersections  with  the  advance 
trajectory,  simultaneously  under  construction,  are  recorded  by  the  program, 
for  they  constitute  the  coordinates  of  the  advance  curve.  When  advance  is 
complete,  subsequent  kinematic-wave  trajectories  either  intersect  the  down- 
stream field  boundary  x = L,  or  trail  off  to  end  spontaneously  at  points  in  the 
x-t  plane  at  which  the  depth  on  the  trajectories  has  fallen,  gradually,  to  zero. 
On  those  trajectories  reaching  x = L,  the  discharge  there  is  recorded,  for  it 
constitutes  the  runoff,  when  the  stream  is  flowing  off  freely,  and  the  rate  of 
ponded-volume  buildup  when  runoff  is  blocked  by  a dike.  The  end  points  of 
those  trajectories  that  trail  off  with  zero  depths  are  also  recorded,  for  they 
constitute  the  coordinates  of  the  recession  curve.  At  the  instant  of  cutoff  of 
the  inflowing  stream,  however,  the  upstream  depth  drops  instantaneously  to 
zero,  in  accord  with  the  normal-depth,  stage-discharge  relation  presumed  in 
force  everywhere  and  for  the  entire  duration  of  the  irrigation.  If  trajectories 
were  constructed  only  for  the  instant  just  before  cutoff  and  the  instant  just 
after,  all  the  detail  of  the  recession  curve  and  some  of  the  runoff  function 
would  be  lost.  Consequently,  the  “instant”  over  which  cutoff  occurs  is  sub- 
divided into  N portions,  and  wave  trajectories  are  computed  for  each  value 
that  the  upstream  depth  takes  on  (at  the  end  of  each  portion  of  the  “instant”) 
as  it  drops  to  zero.  Since  the  instant  is  indeed  infinitesimal  in  size,  all  these 
trajectories  emanate  from  the  single  point  in  the  x-t  plane  (0,  tco). 

When  all  trajectories  have  been  computed,  either  to  their  intersections  with 
the  advance  curve  or  to  the  field  end  or  to  their  trailing  ends  where  y = 0,  as 
described  above,  the  streamflow  calculation  is  over.  If  the  stream  is,  in  fact, 
ponded  behind  an  end  check,  the  infiltration  of  the  still  pool  is  computed 
next.  Finally,  the  program  skips  to  Post-Irrigation  Analysis  (sec.  5.5). 


5.4  SOLMOD  = 5 
Hybrid:  Zero-Inertia  plus 
Kinematic  Wave 


The  great  disparity  in  recessions  computed  by  kinematic-wave  and  zero- 
inertia  theories,  especially  near  the  upper  end  of  the  field,  as  discussed  in 
section  3.16,  led  to  the  development  of  a hybrid  model.  Zero-inertia  theory  is 
used  until  both  advance  is  complete  and  recession  has  started. 


Then,  N kinematic-wave  trajectories  are  constructed  emanating  from  the  node 
points  on  this  time  line  and  utilizing  the  existing  values  of  depth  there  as 
starting  values  (in  case  the  last  depth  profile  computed  by  zero-inertia  theory 
is  highly  irregular,  it  is  smoothed  somewhat  prior  to  passage  to  kinematic- 
wave  theory).  Node  points  are  utilized  in  order  from  left  to  right.  As  in  the 
description  of  the  pure  kinematic-wave  model  (SOLMOD  = 3),  intersections 
of  the  kinematic-wave  trajectories  with  field  end  yield  runoff  or  ponded- 
volume  buildup  rates.  The  trailing  ends  of  those  trajectories  that  fail  to  reach 
field  end  comprise  points  on  the  recession  curve.  When  all  kinematic-wave 
trajectories  have  been  computed,  and  any  still-pool  infiltration  calculated,  the 
program  skips  to  Post-Irrigation  Analysis  (sec.  5.5). 
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5.5  Post-Irrigation 
Analysis 


The  advance  and  recession  curves  computed  at  each  time  step  or  for  each 
kinematic-wave  trajectory  are  printed,  if  IDIAG  > 1,  and  are  used  to  deter- 
mine the  infiltration-opportunity  time  at  n -t-  1 equally  spaced  locations  in  the 
border  (n  is  defined  in  the  sequel).  This  together  with  the  given  infiltration 
function  yields  the  post-irrigation  infiltration  profile  and  the  total  volume  in- 
filtrated. Smooth  profiles  need  only  a small  value  of  n for  proper  definition, 
while  highly  irregular  profiles  require  a iarge  value.  Starting  with  n = 10,  the 
volume  is  computed  and  compared  with  that  stemming  from  a doubled  value 
of  n.  If  the  two  computed  values  of  volume  substantially  agree,  the  larger 
value  of  n is  accepted.  Otherwise  n is  doubled  again.  In  any  case,  n is  not 
allowed  to  exceed  80;  if  it  appears  that  a greater  value  is  necessary,  a warn- 
ing statement  is  printed,  but  the  results  with  n = 80  are  accepted. 

The  ultimate  infiltration  profile  and  advance  and  recession  trajectories  are 
printed  in  dimensionless,  metric,  and  English  units.  If  IDIAG  = 0,  only  11 
equally  spaced  values  are  printed;  with  higher  levels  of  diagnostics,  all  n -i-  1 
computed  values  are  printed. 

This  is  followed  by  some  indicators  of  the  goodness  of  the  computation.  The 
principal  measures  of  such  merit  are  the  degree  of  satisfaction  of  overall 
mass  balance  and  the  smoothness  of  the  computed  surface-water  profiles. 
Dimensionless  values  of  inflow  volume  Vq,  infiltrated  volume  V^,  and  runoff 
volume  V,o  are  printed  as  well  as  the  relative  volume  error  = (Vq  - - VJ/Vq. 

The  surface-water  profiles  were  scanned  during  the  period  of  their  construction 
and  the  largest  value  of  TOOTH  (eq.  4.1)  encountered  is  stored  as  TOOTHM.  At 
the  same  time,  the  largest  relative  vaiue  of  saw  tooth  RTOOTH,  defined  as 

RTOOTH  = T00TH/y3,g  (5.1) 

with 

Yavg  = (Yk  + i + 2yk  + yk_i)/4  (5.2) 

is  stored  as  RTOOTM.  The  smoothness  of  the  computation  is  described  by  the 
following  printed  variables 

ITM  : the  time  index  at  which  the  maximum  value  of  TOOTH  is 

encountered 

KTM  ; the  cell  index  pertaining  to  this  maximum 

TOOTH  : the  maximum 

IRTM  : the  time  index  at  which  the  maximum  relative  value  RTOOTH 

is  encountered 
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KRTM  : the  cell  index  pertaining  to  this  maximum 

RTOOTM  : this  relative  maximum 

The  figures  of  computational  merit  are  followed  by  a synopsis  of  the  results  of 
the  irrigation.  If  IPRZA  = 1 (line  11),  this  commences  with  the  monotonic- 
increasing  infiltration-depth  distribution  derived  from  the  computed  ultimate  in- 
filtration profile,  together  with  the  corresponding  total  infiltrated  volume,  all  in 
dimensionless  units. 

The  remaining  factors  presented  in  the  synopsis  are  given,  where  applicable,  in 
English,  metric,  and  dimensionless  units,  as  well  as  in  percentage  of  the  ap- 
plied depth.  The  variables  presented  are: 


tco 

k 

^FR 

tE 

Vq 

V. 

V,o 

Ymaxu 


Ymax 

^Vinax 


*LQ 

^req 

Zq 

Zavg 


: time  of  cutoff 
: duration  of  advance 
: time  recession  starts  at  upstream  end 
: time  recession  starts  at  downstream  end 
: time  all  surface  water  disappears 
: applied  (inflow)  volume 
: ultimate  computed  infiltrated  volume 
: total  computed  runoff  volume 

: maximum  depth  of  surface  stream  encountered  at  upstream 
end  of  field 

: maximum  depth  of  surface  stream  encountered  at  down- 
stream end  of  field 

: maximum  depth  of  surface  stream  encountered  anywhere 
: location  of  y,^^, 

: minimum  post-irrigation  infiltration  depth 
: maximum  post-irrigation  infiltration  depth 
: average-low-quarter  depth  of  infiltration 
: inputted  required  depth  of  infiltration 
; average  applied  depth  = Vq/L 
: average  infiltrated  depth  = VJL 
: average  depth  of  runoff  = V,o/L 


The  next  three  parameters  comprise  average  depth  z^p  of  deep-percolation 
volume  V.p, 


(5.3) 


with  deep  percolation  based  on  different  definitions  of  required  depth, 
as  follows: 


^Pzreq 

^^Pzmin 


Z 


dPzLQ 


: deep  percolation  defined  as  infiltration  in  excess  of  given  in- 
put value,  z,ep 

: deep  percolation  defined  as  infiltration  in  excess  of 
minimum  infiltration  depth 

: deep  percolation  defined  as  infiltration  in  excess  of  average- 
low-quarter  infiltration  depth 


This  is  followed  by  three  conceptions  of  average  useful  depth  of  infiltration  z^ 


(5.4) 


in  which  is  the  useful  infiltrated  volume  (comprised  of  infiltration  depths 
less  than  or  equal  to  the  required  depth),  namely, 

with  deep  percolation  based  on  inputted  z,eq 

Zu^^.^  : Zu  with  deep  percolation  based  on  minimum  infiltration  depth 

Zu^i^Q  : z„  with  deep  percolation  based  on  average-low-quarter  depth 

of  infiltration 

These  are  followed  by  the  Christiansen  uniformity  coefficient. 


UCc  = 1 


n k=i 


with 


(5.5) 


Zk 


Zk  "I"  ^k+  ^ 
2 


(5.6) 
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the  HSPA  uniformity  coefficient, 


UCh  = 


the  distribution  uniformity, 


DU  = ^ 

^avg 


the  low-quarter  distribution  uniformity. 


DUlq 


and  the  runoff  fraction  in  percent. 


(5.7) 


(5.8) 


(5.9) 


(5.10) 


The  remaining  parameters  are  also  presented  in  terms  of  the  three  interpreta- 
tions of  required  depth— the  inputted  value  z^gq,  the  minimum  in  the  distribu- 
tion z^in,  and  the  average-low-quarter  depth  in  the  distribution  Zlq.  These  com- 
prise the  irrigation  efficiency. 


IE  = 100  — 
Zq 

the  useful  fraction  of  the  infiltrated  volume. 


the  storage  efficiency. 


Zfeq 


(5.11) 


(5.12) 


(5.13) 


in  which  z,eq  stands  for  each  of  the  three  interpretations  in  turn;  the  percent 
of  total  area  adequately  irrigated. 


AAP  = 100 


(5.14) 
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in  which  x^i  is  the  total  length  of  border  receiving  less  than  the  required  depth 
of  infiltration.  This  figure  is  obtained  from  the  monotonic  increasing  equiv- 
alent distribution  of  infiltrated  depths  by  noting  the  x-value  therein  at  which 
the  infiltrated  depth  equals  the  requirement  (interpolating  as  necessary 
between  tabulated  values). 

The  last  parameter  of  the  group  is  the  deficiency  ratio,  DR,  the  average 
deficit  in  the  underirrigated  portion  of  the  field,  expressed  as  a percentage  of 
the  required  depth. 


DR  = 100  ( 1 — — ) 

V y • • 7 * 

^ui  ^req 


(5.15) 


Here,  V^i  is  the  volume  infiltrated  in  the  underirrigated  portions  of  the  border. 
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6.  Equilibrium  (Zero-Inertia)  Model 


The  program  operations  outlined  in  section  5.2  are  detailed  in  this  chapter. 
The  first  nine  sections  provide  the  underlying  theory;  the  rest  sets  forth  the 
program  logic  constituting  the  mathematical  model  embodying  that  theory. 


Part  I:  Theory 


6.1  The  Equilibrium 
Concept:  Subdivision  of 
the  Stream  into  Cells 


In  a useful  approximation  to  reality,  the  forces  on  a slice  of  the  surface 
stream  due  to  the  component  of  weight  acting  down  slope,  the  drag  of 
vegetation  and  soil  surface  acting  upstream,  and  the  pressure  gradient 
resulting  from  nonuniform  depth  can  be  assumed  in  equilibrium.  In  reality, 
any  imbalance  among  these  forces  would  result  in  momentum  changes  in  the 
surface  stream.  In  most  cases,  however,  water  velocities  are  so  low  that  the 
momentum  of  all  portions  of  the  stream  is  nearly  zero,  and  changes  therein 
truly  negligible. 


The  flow  in  the  surface  stream  is  analyzed  by  subdividing  its  length  into  a 
number,  N,  of  cells  and  applying  the  principles  of  conservation  of  mass  and 
equilibrium  of  forces  to  each  of  the  cells.  If  the  shape  of  the  depth  and 
discharge  profiles  were  known,  subdivision  would  be  unnecessary,  and  these 
physical  principles  could  be  used  directly  to  find  the  flow  depths  and 
advance  rate  or  any  other  property  of  the  stream.  The  subdivisions  (cells)  are 
taken  short  enough,  so  that  variation  of  depth  and  discharge  over  the  length 
of  a cell  can  be  assumed  known,  usually  linear.  The  physical  principles  are 
applied  to  these  cells  over  time  steps  short  enough  for  the  shape  of  the  time 
variation  also  to  be  known  (usually,  again,  linear).  The  result  is  a set  of 
mathematical  relationships  among  the  values  of  the  pertinent  variables  at  the 
juncture  sections.  Solution  of  the  equations  yields  the  size  and  shape  of  the 
stream  at  all  times  and,  by  the  way,  the  trajectories  of  the  stream  boundaries 
in  the  x-t  plane. 


36 


Because  the  surface  profiles  are  steep  and  relatively  strongly  curved  when 
the  stream  is  short,  and  vary  more  gradually  as  it  lengthens,  it  is  desirable  to 
have  small  cells  at  small  times  and  larger  cells  at  large  times.  It  is  conven- 
ient to  keep  the  number  of  cells  in  a profile  constant,  from  time  step  to  time 
step  as  pictured  in  figure  6.1.  Consequently,  the  locations  of  the  juncture  sec- 
tions are  moved  over  the  duration  of  a time  step,  as  shown.  In  principle,  a cell 
of  the  surface  stream  is  bounded  by  the  juncture  cross  sections  U and  D on 
its  upstream  and  downstream  sides,  respectively,  by  the  soil  surface,  and  by 
the  moving  free  surface,  so  that  the  cell  is  always  filled  with  water. 

In  figures  6.1  and  6.2,  the  location  of  U at  the  beginning  of  the  time  step  is  Xj, 
at  the  end  of  the  time  step,  XlI  the  downstream  face  is  moved  from  x^  to  Xr. 


Figure  6.2. — Water  cell  with  moving  faces. 
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6.2  Mass  Conservation 
Within  Ceils 


It  is  convenient,  when  expressing  mass  conservation,  to  view  the  cell  volume 
as  consisting  of  both  the  water  above  ground  and  that  which  has  infiltrated 
into  the  soil.  The  cross  sections  “U”  and  “D”  are  simply  extended  below  the 
ground  surface.  Thus,  the  lower  boundary  of  the  cell  lies  below  the  soil  sur- 
face by  an  amount  equal  to  the  depth  of  infiltration  there;  that  is,  by  the 
volume  infiltrated  per  unit  area  of  infiltrating  surface. 

Water  enters  or  leaves  the  cell  through  the  juncture  sections  by  virtue  of  the 
relative  velocity  between  water  and  juncture  section.  The  rate  at  which 
water  actually  enters  the  cell  across  the  moving  upstream  face  U (as  opposed 
to  the  absolute  value  Q,  say  Qj  at  the  beginning  of  the  time  interval)  is 
given  by 


Qr,  = (Vu  - Wu)  Au  - Wu  (6.1) 

Here,  Vy  is  the  water  velocity  at  the  location  occupied  by  U,  as  seen  by  a sta- 
tionary observer,  Wy  is  the  velocity  of  movement  of  U,  Ay  is  the  cross  sec- 
tional area  of  the  surface  stream  at  the  location  of  U,  and  A^^  is  the  volume 
infiltrated  per  unit  length  of  channel  at  U.  Inherent  in  equation  6.1  is  the 
assumption  that  infiltrated  water  moves  solely  in  the  vertical  direction. 

Similarly,  0,^  is  the  rate  at  which  water  flows  out  from  the  cell: 

Qry  = (Vq  ~ Wy)  Ay  — Wy  A^^  (6.2) 

with  the  symbols  defined  as  before;  but  for  the  location  of  the  downstream 
juncture  section  D. 

Mathematical  expression  of  masfe  conservation  of  an  incompressible  liquid  in 
the  given  cell  over  a time  increment  is  given  by 

( ' Q dt  - ( ' Q dt  = ( " A(x,ti)dx  -h  ( ""  A,(x,ti)dx 

•'<1-1  Jxy  JxL  (0-'3) 

’'M  f XM 

A(x,ti_i)dx  - A^(x,ti_i)dx 

Xj  Xj 

The  first  and  second  integrals  on  the  left  represent  the  volumes  of  inflow  and 
outflow,  respectively;  the  first  two  integrals  on  the  right  are  the  volumes  of 
surface  and  infiltrated  water,  respectively,  in  the  cell  at  the  end  of  the  time  in- 
crement; the  last  two  represent  these  same  volumes  at  the  beginning  of  the 
time  increment. 
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There  remains  to  approximate  these  integrals  by  algebraic  expressions  involv- 
ing depth  and  discharge  at  the  node  points  J,  M,  L,  R of  figure  6.1.  The  x - 
integrals  are  expressed  as  follows 

J " A(x,ti)dx  = Al  -I-  0yp  Ar)  (Xr  - Xl)  (6.4) 

J " A,(x,ti)dx  = A,^  -H  </>,p  A,p)  (Xr  - Xl)  (6.5) 

A(x,ti-i)dx  = 6Vjm  = (</>yj  Aj  -I-  (j)y^  Ay)  (Xm  — Xj)  (6.6) 

Az(x,ti_i)dx  = = (0,j  A,j  -I-  <t>,^  A,J  (Xm  - Xj)  (6.7) 

’‘J 

in  which  the  </>  are  shape  factors,  which  depend  upon  the  assumed  variation 
of  the  areas  with  distance.  Usually,  in  most  cells,  a straight-line  variation  of 
the  variables  can  be  assumed,  so  that, 

<^yL  = K = K = ^ 

When  a straight-line  variation  is  considered  inappropriate,  such  as  in  the 
leading  cell  of  the  stream  during  advance,  these  factors  are  changed  accord- 
ingly (see  sec.  6.5). 

The  time  integrals  in  the  left  side  of  eq.  6.3  are  similarly  approximated 
through  a weighting  factor,  d,  but  first,  an  average  face  velocity  w is  expressed 
in  terms  of  face  location  at  the  beginning  and  end  of  the  time  interval.  With  6t 
chosen  small  enough  so  that  the  velocity  of  stream  advance  (or  recession) 
can  be  assumed  essentially  constant  over  the  time  interval,  the  upstream- 
face  velocity  can  be  approximated  by 


(6.9) 


with  a similar  expression  for  Wp.  (One  recalls  that  the  stream  length  is  divided 
into  N cells,  so  that  the  cell  faces  move  in  conjunction  with  the  advancing 
[and  receding]  front  and  rear  stream  boundaries.)  Then  one  can  write 

^ Q,^dt  = I = 0[Ql  6t  - (Al  A,^)  (Xl  - Xj)]  -H  (1  - 0)  (6.10) 

[Qj  8\  - (Aj  -H  A,j)  (Xl  - Xj)] 
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and 


Q,^dt  = O = 0[Qr  6t  - (Ap  + A,^)  (Xp  - Xm)]  + (1  - 0)  (6.11) 

[Qm  — (Am  + (Xp  — Xm)]. 


The  factor  d is  taken  slightly  greater  than  the  value  one-half,  corresponding  to 
straight-line  variation  of  the  integrands  over  time.  This  is  for  reasons  of  com- 
putational stability  as  the  equations  are  applied  successively  over  a long 
sequence  of  time  steps,  rather  than  for  more  accurate  representation  of  the 
average  value  of  the  integrands. 

Equation  6.3  can  thus  be  approximated  by  the  algebraic  equation 


(6.12) 


In  which  is  a residual,  discussed  subsequently. 

In  a border,  with  Q,  A,  and  A^  given  per  unit  width  by  q,  y,  and  z,  respectively, 
the  face  values  of  depth  and  discharge  at  the  node  points  J,  M,  L,  R are  thus 
related  by  equation  6.12,  with 


I = [0qL  + (1  - 0)qj]  5t  - 5Xjl  [d(y^  - zj  -i-  (1  - 0)(yj  -l-  Zj)]  (6.13) 
O = [0qp  -I-  (1  - 0)qM]  6t  - (5Xmr  [0(yR  Zp)  -I-  (1  - 0)(yM  + Zm)]  (6.14) 


= (0y^  yL  + 0yR  yR)  6Xlr 
= (<^ZL  Zl  + </>Zp  Zr)  6Xlr 


(6.15) 


(6.16) 


yj  + ^VM  Vm)  SXjm 


(6.17) 


^^zjM  = Zj  + </>ZM  Zm)  5Xjm 


(6.18) 


in  which 


6Xjl  = Xl  - Xj 


(6.19) 


(6.20) 


5Xmr  — Xp  — Xm 


(6.21) 


(6.22) 
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With  nodal  locations  J,  M,  L,  R,  in  the  x - t plane  given  and  advance  up  to 
time  ti_i  known,  the  infiltration  times  tj,  tm,  tl,  tr  are  known.  In  a uniform 
border,  z depends  virtually  on  r alone,  (is  virtually  independent  of  flow  depth), 
so  that  Zj,  Zm,  Zl,  Zr  are  also  known.  Assuming  that  the  flow  has  been  solved 
up  to  time  ti_i,  depth  and  discharge  at  J and  M are  also  known.  The  equation 
6.12  contains  the  four  unknowns  yL,  Ql,  yR,  Qr.  R,.  represents  a residual  that 
results  from  incorrect  values  of  these  variables;  correct  values  reduce  R^  to 
zero.  An  additional  equation  relating  the  same  four  variables  is  provided  by 
the  statement  of  equilibrium  amongst  all  forces  acting  on  the  cell  at  time  V 

6.3  Equilibrium  of  Figure  6.3  shows,  schematically,  the  longitudinal  components  of  forces  act- 

Forces  Acting  on  a Cell  j^g  on  the  surface  water  contained  in  the  cell  at  time  t|.  The  forces  of  hydro- 
static pressure  exerted  by  the  water  in  the  neighboring  cells  upstream  and 
downstream  from  the  given  cell,  divided  by  the  unit  weight  of  water,  are  given 
by  Pl  and  Pr,  respectively.  W is  the  component  of  the  weight  acting  in  the 
downstream  direction,  and  D is  the  drag  of  soil  surface  and  vegetation  acting 
upstream,  each  also  divided  by  the  unit  weight  of  water.  With  P the  resultant 
P L - Pr  of  the  pressure  forces,  the  statement  of  equilibrium  is  simply 

P -f  W - D = 0 (6.23) 


In  a unit  width  of  border,  the  resultant  of  the  pressure  forces  is  given  by 


(6.24) 


Figure  6.3.  — Equilibrium  of  forces  on  suface- water  cell. 
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The  weight  component  W in  a unit  width  of  border  is 


W = Yr)  (b,  - bn)  (6.25) 

in  which  bL  and  bp  are  the  bed  elevations  at  Xl  and  Xp,  respectivelY- 

The  drag  d per  unit  length  of  stream  (divided  bY  the  unit  weight  7 of  water), 
called  here  the  drag  intensitY,  depends  upon  the  velocitY  and  depth  of  flow. 
The  relationship  is  put  in  concrete  terms  through  the  use  of  a dimensionless 
drag  coefficient  Cq  defined,  in  effect,  bY  the  following  equation 

VIVI 

7d  = w Cd  e (6.26) 

in  which  w is  the  wetted  perimeter  of  the  flow,  unitY  in  a unit  width  of 
border,  and  q is  the  mass  densitY  of  water.  The  absolute-value  signs  are  put 
on  the  velocitY  V to  insure  that  the  drag  is  alwaYS  in  a direction  opposite  that 
of  the  flow  velocitY-  The  drag  coefficient  can  be  expressed  in  terms  of  the 
well  known  ChezY  as  follows.  First  a definition  is  given  to  the  so-called 
friction  slope  S,  through  which  Ch  is  usuallY  defined.  S,  is  the  drag  per  unit 
length  divided  bY  the  weight  of  the  stream  per  unit  length,  the  latter  being  the 
product  of  unit  weight  and  cross-sectional  area.  In  SYmbols, 

S,  = (6.27) 

A 

The  ChezY  Ch  is  usuallY  defined  in  terms  of  the  friction  slope,  bY  the 
expression 


Sf 


V|V| 

Ch^R 


in  which  R is  the  hYdraulic  radius 


w 

It  follows  that  Cd  and  the  ChezY  Ch  are  related  as  follows: 


in  which  g is  the  ratio  of  unit  weight  to  mass  densitY- 


(6.28) 


(6.29) 


(6.30) 
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A variety  of  empirical  expressions  have  been  used  to  predict  the  value  of  Ch, 
which,  in  fact,  depends  upon  the  geometry  of  the  roughness  elements  or 
vegetation,  the  cross-sectional  dimensions  of  the  stream,  and  the  Reynolds 
number  of  the  flow.  The  program  choices  are  discussed  in  section  3.4.  The 
most  common  formula  in  use  in  the  United  States  expresses  the  Chezy  in 
terms  of  the  Manning  n 


C,. 

Ch  = R"'®  (6.31) 

in  which  is  a units  coefficient  (C^  = 1.0  m^'^/s  in  the  metric  system; 

Cu  = 1.486  fV’^lsec  in  the  English  system). 

In  any  event,  the  total  drag  on  the  water  in  the  cell  (divided  by  the  unit 
weight)  is 


D = 


dx  = d 6Xlr 


(6.32) 


in  which  d is  defined  as  the  average  drag  intensity.  It  is  not  a trivial  matter  to 
express  the  average  drag  intensity  in  terms  of  its  values  at  L and  R.  With 


Q|Q|  _ V|Vlw 
ACh^R  " Ch^ 


(6.33) 


small  values  of  flow  depth  make  the  denominator  very  small,  with  the  result, 
of  possibly  very  large  values  of  d and  very  nonuniform  variation  of  d over  the 
cell  length.  Two  variants  are  proposed  for  estimating  d in  a border,  now  com- 
puted per  unit  width.  The  first,  for  a cell  in  which  velocity  varies  gradually,  as 
in  the  lead  cell  during  advance 


VIVI 

d = 

(y) 

Here  V is  the  average  velocity  in  the  cell,  approximated  by 

V = <f)y^  Vl  -H  <^Vr  Vr 


(6.34) 


(6.35) 


with  most  commonly 


1 

</>vl  = <Avr  = y 


(6.36) 
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The  function  is  evaluated  at  an  average  value  of  depth 


y = Vl  + Yr  (6.37) 

in  which  the  0d  are  shape  factors  pertinent  to  the  calculation  of  drag.  In  in- 
terior cells,  with  neariy  straight-line  variation  of  y over  a cell  length 

4>d  = = — (6.38) 

^yl  Yr  2 


The  second  variant,  especially  suited  to  cells  in  which  discharge  varies 
gradually  over  a cell  length,  such  as  in  the  downstream-most  cell  during  run 
off,  or,  for  that  matter,  in  the  interior  of  the  flow,  is 


with 


3 = ^ 

KD(y) 


Q — <^qL  dL  + 4>qf^  ^R 


Here,  commonly. 


1 

0qL  = <AqF|  = — 


(6.39) 


(6.40) 


(6.41) 


The  term  Kp  (y)  represents  the  function  of  y,  evaluated  at  an  average 

value  of  depth,  as  given  by  equation  6.37. 

The  companion  equation  to  (12)  expressing  equilibrium  of  forces  is  then 


R,,=  P-hW-D  = 0 (6.42) 

Like  Rc,  Rm  is  a function  of  the  unknown  variable  values  qL,  yL,  qR,  yp.  It  is  evi- 
dent that  between  them,  the  two  equations,  12  and  42,  contain  four  variables 
and  are  in  themselves  insufficient  for  solution.  Application  of  these  equations 
to  N cells  yields  2N  equations  and  2N  -i-  2 unknowns.  Additional  information 
for  the  boundary  cells  at  front  and  rear  closes  the  system  by  supplying  two 
more  values  of  unknowns. 


6.4  Boundary  Conditions  At  the  upstream  boundary  of  the  stream,  the  discharge  is  presumed  known 

for  all  time  of  interest.  The  inflow  discharge  hydrograph  constitutes  the  up- 
stream boundary  condition.  Most  commonly,  the  inflow  rate  is  constant  until 
cut  off  and  is  zero  thenceforth. 

The  downstream  boundary  condition  depends  upon  the  progress  of  the  irriga- 
tion and  upon  the  assumed  nature  of  the  downstream  end  of  the  border  or 
furrow;  that  is,  whether  it  is  free  flowing  or  diked  by  a check. 

During  advance,  the  depth  and  discharge  at  the  wave  front  are  both  zero. 
However  the  x-position  of  the  wave  front  is  not  known.  Yet  the  preceding  set 
of  algebraic  equations  is  simpler  to  solve  if  the  locations  Xl,  Xr  are  known  a 
priori  for  all  cells.  So,  only  Xr  of  the  most  downstream  cell  is  carried  in  the 
equations  as  an  unknown.  The  interior  values  are  set  by  estimating  the  incre- 
ment of  advance  occurring  in  the  given  time  step  (by  extrapolating  the  advance 
in  the  previous  time  step)  and  dividing  the  estimated  stream  length  into  cells. 
If  the  cell  division  is  not  too  fine  at  the  front,  and  the  time  steps  are  not  so 
great  as  to  cause  significant  error  between  estimated  and  ultimately  deter- 
mined stream  advance,  the  potential  problem  of  negative  cell  lengths  does 
not  arise.  If,  indeed,  a negative  cell  length  is  generated,  the  stream  is  redivided 
into  new  cells. 

The  boundaries  of  the  stream  in  recession  are  set  as  follows:  At  a stream 
boundary  at  which  recession  is  imminent,  the  discharge  is  zero.  Recession 
from  the  boundary  location  during  a given  time  step  is  signalled  by  the  com- 
putation of  a negative  depth  there.  When  this  happens,  the  computational 
stream  boundary  is  moved  into  the  interior  of  the  stream  until  a location  is 
found  at  which  the  depth  exceeds  some  minimum  small  value  y^ec  (perhaps 
1 pet.  of  the  normal  depth  for  the  initial  inflowing  discharge).  This  constitutes 
the  new  location  Xlb  (or  Xrb,  if  the  recession  occurs  from  the  downstream  end) 
of  the  computational  stream  boundary.  The  discharge  there  is  set  to  zero, 
with  the  implication  that  all  water  outward  from  that  point  is  stagnant  and 
simply  seeps  vertically  into  the  soil. 

Prior  to  any  recession  from  the  downstream,  or  front,  end  of  the  surface 
stream,  one  of  three  conditions  prevails  there.  Either  the  stream  is  advancing, 
as  discussed  previously,  or  it  is  running  off  the  field  into  a drainage  ditch  or 
it  is  ponded  behind  a dike.  Other  conditions  are  possible,  for  example,  flow 
over  a dike  of  insufficient  height  to  contain  the  stream,  but  these  have  not 
been  programmed.  With  a dike  in  place  at  field  end,  the  discharge  there  is 
simply  zero. 

In  the  event  of  free  overfall  into  a drainage  ditch,  the  depth  of  flow  is, 
perhaps  startlingly,  zero.  This  conclusion,  supported  in  Strelkoff  and 
Katopodes  (9),  stems  directly  from  the  assumption  of  zero  inertia  in  the  flow. 
In  the  usual  hydraulic  treatment  of  subcritical  flow  approaching  a free  over- 
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fall,  it  is  noted  that  the  water  accelerates  in  the  vicinity  of  the  brink,  until  the 
depth  drops  to  critical,  at  which  section  the  specific  energy  of  the  flow  is  a 
nninimum,  and  the  velocity  head  is  a given  fraction  of  the  flow  depth,  the  frac- 
tion depending  on  the  cross  section  (one-half,  with  a rectangular  section).  In 
section  6.8,  the  assumption  of  zero  inertia  without  zero  weight  is  tantamount 
to  assuming  a ratio  of  mass  density  to  unit  weight  of  zero,  or  an  infinite  value 
of  g.  Thus,  all  velocity  heads  are  zero,  relative  to  flow  depths;  all  flows  are 
necessarily  infinitely  subcritical.  The  ratio  of  critical  depth,  to  normal  depth, 
or  any  other  depth  in  the  flow  is  thus  zero. 

The  result  that  velocity  at  the  brink  is  infinite  is  dealt  with  from  several  points 
of  view  by  Strelkoff  and  Katopodes  (9).  From  the  theoretical  standpoint,  this 
result,  obviously  In  gross  error,  stems  from  strict  adherence  to  the 
equilibrium  assumption,  even  in  the  vicinity  of  the  brink,  where,  in  fact,  the 
forces  are  not  in  equilibrium.  To  a lesser  degree,  the  same  problem  arises 
with  the  use  of  the  Saint-Venant  equations  or  the  equation  of  gradually  varied 
steady  flow  near  a brink.  The  latter  formula. 


dy  _ So  - S( 
dx  ~ 1 - F2 


(6.43) 


in  which  F is  the  local  Froude  number  of  the  flow 

V 

F = - y . 

Vg  A/B 


(6.44) 


dv 

with  B the  top  width,  yields  the  result  — ^ = oo  at  the  critical  section,  where 

dx 


F = 1.  Now  observations  of  the  phenomena  in  real  life  show  that 


dy 

dx 


is  not 


infinite  at  all,  but  rather  small,  well  less  than  unity.  The  grossly  incorrect 
results  of  equation  6.43  stem  from  the  strict  adherence  to  the  assumption  of 
hydrostatic  pressure  distribution,  even  when  it  does  not  apply,  such  as  close 
to  the  brink  where  the  streamlines  are  distinctly  curved. 


If  it  were  necessary  to  model  the  flow  in  the  vicinity  of  the  brink  accurately,  a 
more  accurate  theory  would  have  to  be  invoked.  In  general  practice,  it  has 
become  customary  to  ignore  the  errors  incurred  by  use  of  equation  6.43  or  its 
counterpart  in  unsteady  flow,  the  Saint-Venant  equation,  because  in  com- 
parison to  the  long  reaches  of  open  channel  usually  considered,  the  length  of 
channel  over  which  significant  errors  arise  is  actually  very  small.  Similarly,  as 
shown  by  numerical  comparisons  in  the  foregoing  reference,  provided  the 
Froude  number  everywhere  in  the  flow  is  in  fact  small,  say,  less  than  about 
0.3,  the  longitudinal  extent  of  the  region  with  abnormally  high  velocity  values— 
computed  in  the  assumption  of  force  equilibrium— is  also  very  small,  relative 
to  the  overall  stream  length.  Thus,  from  a practical  standpoint,  one  can  sim- 
ply ignore  the  very  high  velocity  values  that  are  computed  for  this  very  small 
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fraction  of  the  stream  length.  From  the  standpoint  of  the  present  computa- 
tional technique,  it  is  not  even  velocity,  but  discharge,  which  constitutes  the 
primary  variable,  and  the  variation  of  discharge  with  distance  is  well  behaved, 
even  near  the  brink.  Only  the  computation  of  average  drag  intensity  (eq.  6.32) 
in  the  cell  at  the  brink  requires  special  care.  Clearly,  right  at  the  brink,  drag 
intensity  is  infinitely  great  (eq.  6.33)  and  is  not  at  all  representative  of  condi- 
tions within  the  cell;  consequently,  equation  6.39  is  used. 


A similar  situation  arises  in  the  leading  cell  during  advance.  While  the  veloc- 
ity variation  is  well  behaved  near  the  front,  the  zero  depth  there  again  results 
in  an  infinite  drag  intensity  again  not  representative  of  average  conditions  in 
the  cell,  and  equation  6.34  is  appropriate. 

The  matter  of  expressing  cell  averages  in  terms  of  end-point  values  is,  in 
general,  not  trivial.  The  question  is  dealt  with  in  some  detail  in  the  next 
section. 


6.5  Cell-Profile  Shape 
Factors:  The  0 
Weighting  Factors 


The  weighting  factors  are  the  key  to  a precise  solution  of  the  governing  equa- 
tions for  the  water  in  the  cells,  and  hence  of  the  entire  stream  profile.  The 
algebraic  cell  equations  are  precise  equivalents  of  the  integral  formulas 
(which  are  exactly  correct)  if  the  weighting  factors  are  chosen  just  right. 
Usually,  the  smaller  the  cell  length,  the  less  troublesome  is  selection  of 
appropriate  values  of  the  weighting  factors.  Economical  calculation,  however, 
suggests  use  of  large  cells,  and  consequently  less  obvious  values  of  the  (t>. 
Poorly  chosen  (j)  can  result  in  the  development  of  saw-tooth  profiles;  this  can 
be  troublesome,  because  of  the  potential  for  calculating  negative  depths  in 
the  interior  of  the  stream.  This  would  not  represent  a true  instability  of  calcu- 
lations, because  small  errors  are  not  amplified  without  bound  as  the  com- 
putations proceed  over  the  succession  of  time  steps.  The  end  result,  how- 
ever, can  be  the  same:  an  aborted  computation. 


As  suggested  in  sections  6.2  and  6.3,  the  variation  of  depth  and  discharge  in 
the  interior  of  the  stream  is  sufficiently  gradual  and  close  to  linear  that  the 
trapezoidal  rule  of  numerical  integration  can  be  readily  applied  even  when  the 
number  of  cells  is  small  and  the  length  of  each,  consequently,  large.  Then, 
the  (j)  weighting  factors  all  have  the  value  one-half  as  in  equations  6.8,  6.36, 
6.38,  and  6.41.  Near  the  boundaries  of  the  surface  stream,  however,  variables 
can  change  rapidly  and  curvilinearly.  This  is  clearly  true  near  the  front  of  an 
advancing  stream  or  near  a free  overfall.  It  is  true,  also,  at  an  end  check  for  a 
short  period  after  the  stream  arrives.  Finally,  conditions  near  the  trailing  edge 
of  a stream  in  recession  on  a steep  slope  may  also  require  special  treatment. 

The  tip  cell  is  investigated  first.  In  the  advancing  tip,  velocity  variation  is 
smooth  and  gradual,  while  discharge  is  not.  In  the  free  overfall,  discharge 
variation  is  smooth  and  gradual,  while  velocity  is  not.  In  both,  the  variation  in 
depth  is  highly  curvilinear.  The  depth  variation  in  the  tip  cell  can  be  approxi- 
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mated  by  assuming  it  identical  to  that  of  a corresponding  steady  flow.  In  the 
advancing  stream,  the  corresponding  steady  flow  is  obtained  by  postulating 
that  the  stream  advances  without  changing  shape,  that  is,  that  all  cross  sec- 
tions move  at  the  same  average  velocity  (that  is,  averaged  over  the  cross- 
sectional  area)VA,  equal  to  the  velocity  of  advance.  Indeed,  the  average  fluid 
velocity  in  the  surface  stream  changes  but  slowly  with  both  time  and 
distance.  The  details  of  flow  at  the  very  front  of  the  wave  are  obscure,  but  it 
is  a reasonable  postulate  that  the  streamlines  intersect  the  wave  front,  as 
shown  in  figure  6.4.  There  is  a velocity  gradient  in  the  surface  stream,  with 
velocity  increasing  with  distance  from  the  bed,  because  of  bed  resistance 
and  infiltration.  This  concept  agrees  with  the  rolling  advance  of  the  front  evi- 
dent in  the  field,  as  opposed  to  bulk  sliding,  which  might  at  first  seem  to  be 
the  consequence  of  the  flow  pattern  postulated  in  figure  6.4.  In  any  event, 
equation  6.33,  with  the  Chezy  Ch  defined  as  in  the  bulk  of  the  wave®,  is 
assumed  to  apply  all  the  way  to  the  stream  front.  The  depth  profile  near  the 
front  is  then  obtained  as  follows. 


Equation  6.23  written  for  a very  thin,  unit-width  slice  of  stream  moving 
downstream  at  the  velocity  is  given  by 


-I- 


Vl  + Vr 
2 


So  6x 


(6.45) 


in  which  is  the  Chezy  Co  for  an  average  depth  in  the  slice.  As  the  slice 
thickness  is  decreased,  the  average  value  of  depth  therein  approaches  the 
depth  y at  midsection.  Division  of  equation  6.45  by  the  volume  of  the  slice 
and  passage  to  the  limit  as  6x  — 0 yields  the  differential  equation 


d 


Figure  6.4. — Assumed  flow  pattern  in  stream  front. 


. .with  the  exception  of  those  formulas  like  the  Sayre- 
Albertson  (6),  which  fail  below  a certain  positive  value  of  y. . . 


(6.46) 


dx 


The  behavior  y(x)  is  readily  obtained  from  equation  6.46  by  numerical  integra- 
tion. To  start  the  process,  it  is  noted  that  for  x — x^,  the  right  side  increases 
without  bound  as  y — 0,  while  So  remains  fixed.  Thus,  very  close  to  the 
stream  front,  Sq  is  negligible,  and  the  differential  equation  becomes 
separable  with  a solution  that  depends  upon  the  assumed  variation  of  C^,  with 
y.  If  the  Manning  formula  (eq.  6.31)  is  used,  the  result  is 


(6.47) 


In  a horizontal  border  this  expression  is  valid  arbitrarily  far  back  from  the 
front  as  long  as  V is  reasonably  close  to  V^.  In  a sloping  border,  equation  6.47 


step  in  a numerical,  stepwise  integration  of  equation  6.46,  say  by  the 
modified  Euler  method. 

The  shape  factor  r^  (4)^^  for  the  tip  cell)  can  be  estimated  by  integrating  y(x) 
over  the  length  of  the  cell,  then  dividing  by  cell  length  to  get  the  average 
depth  yLR.  The  ratio  of  average  depth  to  depth  yN  at  the  left  boundary  of  the 
tip  cell  is  the  shape  factor,  that  is. 


(6.48) 


since  yN^,  is  zero.  Evidently,  the  shape  factor  will  vary  with  the  length  of  the 
tip  cell.  For  a very  short  cell,  in  view  of  the  3/7  power  of  the  variation,  the 
shape  factor  is  simply. 


With  increasing  length,  this  factor  increases.  A unique,  generalized  solution 
for  the  depth  profile  near  the  front  of  a stream  advancing  with  constant 
velocity  on  a sloping,  impervious  bed  can  be  developed.  The  variation  of  y* 

( = y/y^)  is  found  as  a function  of  x^*  - x*  ( =[Xa  - x]So/yn),  in  which  y^  is 
the  normal  depth  for  the  given  velocity  of  advance  and  given  border  slope  and 
roughness.  (See  sec.  6.8  for  further  discussion  of  dimensionless  parameters.) 
Consequently  the  shape  factor  will  also  be  a unique  function  of  dimension- 
less tip-cell  length  6xn*.  This  relationship  was  computed  and  tabulated,  and 
an  empirical  expression  was  devised  to  approximate  the  relation,  namely. 


Tg  = 1.0  - 0.275 


(6.50) 
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When  needed,  as  in  the  passage  from  the  first  time  step  of  solution  to  sub- 
sequent time  steps  (sec.  6.10),  depth  values  can  be  obtained  at  points  in  the 
interior  of  the  end  cell  by  postulating  a power-law  variation  of  depth  with 
distance  from  the  brink;  the  power  ^ used  is  that  necessary  to  yield  the  given 
shape  factor. 


1 


with,  then, 

Yn  'Xa  - Xn' 


(6.51) 


(6.52) 


Here,  is  the  location  of  the  left  end  of  the  tip  cell. 

A similar  procedure  is  invoked  to  estimate  the  shape  factor  in  the  tip  cell 
under  run-off  conditions.  Only,  now,  the  downstream  boundary  of  the  cell  is 
fixed,  and  the  discharge  (over  an  impenetrable  surface),  rather  than  the 
velocity,  is  constant,  at  the  run-off  rate  q,o.  The  appropriate  differential  equa- 
tion, replacing  equation  6.46,  is 


dx  ° 0^2  y" 


(6.53) 


The  resuiting  power-law  solution  at  vaiues  of  L - x so  small  that  So  is 
negligible  is 


y = 


(L  - x) 


(6.54) 


in  which  L is  the  field  length  (x-coordinate  of  field  end);  the  Manning  formula 
for  resistance  was  again  employed,  to  obtain  this  result.  Equation  6.54,  thus, 
replaces  equation  6.47.  The  shape  factor,  r„  for  the  run-off  case. 


r„  = 


Vlr 

Yn 


is  readily  found  for  a very  short  cell  length,  namely 

13 
r„  = 


16 


(6.55) 


(6.56) 
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Numerical  integration  of  equation  6.54,  integration  of  the  results  to  yield  Ylr 
and  hence  r„  as  a function  of  L*  - x*^,  and  fitting  an  empirical  expression  to 
the  last  yields  the  approximate  relation  for  longer  cells,  namely 

r„  = 1 - 0.12  (6.57) 

in  place  of  equation  6.50.  Again,  the  exponent  in  an  approximately  fitting 
power-law  relation  for  depth. 


a = — - 1 (6.58) 

Profiles  in  the  vicinity  of  the  blocked  end  of  a border  can  also  exhibit  highly 
irregular  behavior,  as  illustrated  by  the  depth  profile  of  figure  6.5.  It  is  possi- 
ble to  make  an  estimate  of  the  location  of  the  reentrant  corner  by  computing 
the  volume  that  would  have  run  off  were  the  dike  absent  and  then  assuming 
that  this  volume  forms  a stagnant  pool  at  the  downstream  end.  This  pro- 
cedure provides  an  approximate  profile  shape,  and  so  can  lead  to  weighting 
factors.  The  complexity  of  the  approach  suggests,  however,  simply  to 
increase  the  number  of  cells  until  the  error  from  this  source  is  reduced  to 
tolerable  levels.  Indeed,  the  program  BRDRFLW  does  this. 


X 


N 


Figure  6.5.  — Interaction  of  live  stream  with  pond. 


9L*  = LSo/Vn- 
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During  advance,  a shape  factor  different  from  one-half  is  also  needed  for  the 
z-profile  in  the  lead  cell,  say  at  time  X,.  The  shape  is  determined  by  assuming 
that  over  a short  time  period  prior  to  t|,  the  stream  had  advanced  at  constant 
velocity  Va-  At  the  instant  t|  then,  any  point  x behind  the  stream-front  location 
Xa  would  have  been  wetted  for  a period  of  time 


T = 


(6.59) 


If  the  Kostiakov  infiltration  formula,  z = kr®,  is  assumed, 

rXA  (ti)  - X V 


z(x,ti)  = k (- 


Va  (ti)  ' 

It  follows  that  the  z-profile  in  the  lead  cell  is  given  by 

Z / Xa  - X \a 


/ - X Y 

' 6Xn  f 


and  the  shape  factor  for  the  lead  cell  is 


(6.60) 


(6.61) 


0z,  = Ta  = 


■ 1 + a 


(6.62) 


In  the  lead  cell  during  advance,  z^  + i = 0.  Formula  6.62  is  valid  for  long  cells, 
as  long  as  the  advance  velocity  has  not  changed  very  much  over  the  period  of 
time  necessary  to  advance  a distance  equal  to  the  length  of  the  cell.  Of 
course,  an  infiltration  equation  differing  from  the  monomial  power  law  above 
can  lead  to  a different  shape  factor. 


Because  the  curvature  of  the  infiltrated-depth  profile  near  its  downstream  end 
persists  for  some  distance  upstream  from  an  advancing  front,  and  also  for 
some  time  in  the  vicinity  of  field  end  after  advance  has  been  completed,  a 
formula  for  the  shape  factors  (j)^^  and  was  devised  for  general  use.  This 
formula,  just  for  appropriate  shape  factors,  is  based  on  the  assumption  that 
the  shape  of  the  profile  in  any  cell  can  be  derived  from  a monomial  power 
law  in  distance  back  from  a fictitious  (or  real)  advancing  front.  The  location  of 
this  front  is  determined  by  given  depths  of  infiltration  Zl  and  Zr  on  the  left 
and  right  sides  of  the  cell,  respectively,  the  cell  length,  6x  = Xr  - Xl,  and  the 
assumption  that  the  power  is  a.  With 


a = 


(6.63) 
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it  follows  that  the  distance  between  Xl  and  the  (real  or  fictitious)  profile 
front  where  the  infiltration  depth  drops  to  zero  is  given  by 

5x 

1 (6.64) 

1 — ff 

Then  with  the  infiltrated  volume  for  the  cell,  given  by  both  the  equations 

6V,  = = (4,  z,  + z„)  SX  (6.65) 

1+a  1+a  ” 

(in  which  - 6x)  and  the  definition 

<^ZR  = 1 - K 

it  follows  that 


0zi  = 

These  general  formulas  show  that  (j)^^  and  both  approach  1/2  as  Zr  — z^_, 
and  that  <t>^^  approaches  1/(1  + a)  (eg.  6.62),  as  Zr  - 0. 


(6.66) 


(6.67) 


6.6  Solution  of  the  Set 
of  Nonlinear  Algebraic 
Equations  by  the 
Newton-Raphson 
Technique:  Successive 
and  Single  Linear 
Approximation 


As  indicated  in  section  6.3,  application  of  mass  conservation  and  force  equili- 
brium to  the  N cells  making  up  the  total  stream  length  yields  2N  equations  of 
the  type  6.12  and  6.42;  namely. 


Rck  = Ik  - O,  - 


Vk 


+ ^X'knownk 


^^^knowrik  “ ® 


(6.68) 


and 


Rm,  = Pk  + W,  - D,  = 0 (6.69) 

with  k = 1,  2,  3,  . . .,  N identifying  the  cell.  These  contain  the  2N  -i-  2 
variables,  yk,  dk,  k = 1,2,  . . .,  N -i-  1.  The  values  yknown,  hknown  are  presumed 
known  and  represent  the  solution  at  the  previous  time  step.  From  the 
previous  discussion  on  boundary  conditions,  Pi  is  known. 


[ q,  = Pin  t < too  ] 

for 

[ Pi  = 0 t > U ] 


(6.70) 
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In  the  simplest  case,  is  the  constant  inflow  until  cut  off.  During 
advance, 


Qn  + 1 — Yn  + I — 0 


(6.71) 


but  x^.,  and  hence  6xn  is  not  known.  The  remaining  cell  lengths  6Xk,  k = 1,  2, 

. . .,  N - 1,  are  known  from  extrapolation  of  the  advance  over  the  previous 
time  step.  With  uniform  division  of  stream  length  into  N parts  (except  for  the 
tip  cell),  and  uniform  6t, 


6Xk 


+ (XAi_,  - Xa^_,)  - X, 


N 


k = 1,2, 


N - 1 (6.72) 


the  upstream  computational  boundary  Xlb  = 0.  provided  recession  has  not 
begun. 


During  run  off. 


Vn.i  = 0 (6.73) 

and  Pn  + i is  unknown,  while  with  ponding  behind  a dike, 

Qn.i  = 0 (6.74) 


and  Yn  + i is  unknown. 

The  nonlinear  nature  of  the  force  equation  makes  it  necessary  to  solve  the 
set  of  equations  and  boundary  conditions  iteratively,  unless  a locally  linear- 
ized solution  is  considered  adequate.  In  principle,  as  long  as  the  change  in 
value  of  a variable  over  a time  step  is  small,  relative  to  the  value  of  the 
variable  itself,  the  locally  linearized  solution  is  accurate  and  avoids  the 
repeated  calculations  of  an  iterative  solution.  Thus,  using  a small  enough 
value  of  time  step  will  usually  be  adequate  to  insure  the  accuracy  of  the 
linearized  solution.  On  the  other  hand,  in  those  regions,  such  as  near  the 
trailing  edge  in  recession  where  y and  q are  both  very  small,  there  is  some 
question  as  to  how  small  the  changes  thereof  have  to  be  to  insure  accuracy 
of  the  linear  solution.  Further,  while  the  nonlinear  solution  requires  iterative 
procedures,  the  time  steps  can  be  larger  than  in  the  linearized  scheme.  Until 
the  user  has  gained  familiarity  with  the  behavior  of  the  mathematical  model, 
it  is  recommended  that  the  iterative,  nonlinear  solution  be  used.  For  only 
then  is  the  user  assured  that  the  equations  of  mass  conservation  and  force 
equilibrium  have  been  satisfied.  In  any  event,  the  linear  solution  comprises, 
simply,  the  first  step  of  the  iterative  solution. 
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It  is  desired,  then,  to  find  those  values  of  Vk,  q,,  as  satisfy  the  inequalities 


RcJ  ^ Rc 

^max 

RmJ  Rm 

"'k'  ^max 


k = 1,  2,  ...,  N 


(6.75) 


in  which  R^  and  R„  are  small  quantities  on  the  order  of  0.000001  times 
the  size  of  the  individual  members  of  the  right-hand  sides  of  equations  6.68 
and  6.69.  Substitution  of  first  guesses,  say  yknownk.  Qknownk.  ^or  the  unknown 
values  yields  residuals  Rcj^,  Rn^^^  generally  greater  than  the  maximum.  So,  the 
guesses  must  be  adjusted  by  an  amount  as  wiil  reduce  the  residual.  With  j 
the  iteration  index,  the  quantities 


^k  = 5yk  = - y'k 

V = 6q,  = qi  + ^ - q'k 


(6.76) 


are  sought,  in  which  the  6 symbol  refers  to  the  change  from  one  iteration  to 
the  next  on  a given  time  line.  They  are  found  by  solution  of  the  simultaneous 
set  of  linear  equations 


dR 


6Ro,  = 


Ck 


dR 


5Rmk  - 


dVK 

3Rmk 

5yk 


^k  + 


Ck 


aR 


^k  + 


9yk+i 

^Rmk 


’k  + 1 


-t- 


Ck 


aR, 


3qk 

aR 


^k  + 


Ck 


ay 


k + 1 


^k+1  + 


rHk 


^*^k  + 1 

aR 


^k+1  = 


■Ck 


aqk 


Vk  + 


rrik 


aq 


Vk+^ 


= -R 


k + 1 


nik 


(6.77) 


The  partial  derivatives  are  determined  for  the  approximation  to  yk  and  Pk 
by  differentiating  equations  6.68  and  6.69  and  substituting  the  values  of  yk 
and  Pk.  The  set  of  equations  6.77  is  linear  and  is  solved  by  the  methods  of  the 
next  section.  For  k = N,  during  advance,  since  y^  + i and  Pn  + i are  both  known, 
i7n  + i is  defined  as  the  correction  to  the  length  of  the  tip  cell 


= xj^+^  - x^  (during  advance)  (6.78) 

and  the  coefficients  of  rj^  + i in  equations  6.77  are,  correspondingly,  dRJdx^^ 
and  dRJdXf^.  These  values  occupy  the  slots  in  the  coefficient  matrix  that 
would  normally  contain  SRc^/^Pn  + i and  ^R^J^qN  + ^,  as  the  latter  are  of  no 
interest. 


With  yL,  Pl  representing  yk,  Pk,  respectively,  and  yR,  Pr,  yk+i,  Pk+i,  respectively, 
the  Newton-Raphson  coefficients  are,  for  each  k, 
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— d 6XjL  — 5XlR 


(6.79) 


3Rc 

dVL 


3Rc 

— ^ = 0 5t  (6.80) 

aqi 

5Rc 

= d 5Xmr  — 0yp  axLR  (6.81) 

3yR 

except  for  k = N in  the  event  that  advance  halted  during  the  current  time 
step.  In  that  case, 


aRc 

= - axLR  (6.82) 

9yR 

because  the  movement  of  the  right  boundary  does  not  “scoop  in”  any 
volume.  Finally, 


aRc 


- a at 


(6.83) 


During  advance. 


axR  SXlr 

all  these  partial  derivatives  stem  from  equation  6.68. 
From  equation  6.69  the  force  equation. 


(6.84) 


aRm  / 1 aCf,  1 \ 

= yL  + K (bL  - bR)  + 2 0n  , D ( ^ — ^ + - ) (6.85) 

dVL  ' Ch  ay  y ' 


Here,  the  values  D and  y are  connputed  as  per  equations  6.32  and  6.34  or 
6.39.  The  form  of  the  derivative  dCJdy  depends  on  the  empirical  expression 
used  to  evaluate  (see  sec.  3.4).  With  the  Manning  formula  in  force. 


dC,  ^ 

ay  6 y 


(6.86) 


in  which  Cn  is  the  Chezy  Ch  evaluated  at  the  average  value  of  y.  Further, 


aqi 


D 


(6.87) 
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Some  problems  were  experienced  with  convergence  of  the  iterative  scheme 
(or  development  of  saw-tooth  profiles  near  the  receding  trailing  edge  with  the 
linear  scheme)  when  this  expression  was  used  for  the  coefficient  of  SQl  (and, 
with  replaced  by  for  the  coefficient  of  Sqp).  These  problems  stemmed 
from  the  fact  that  near  the  trailing  edge  in  recession,  q values  are  very  small. 
Equations  6.77,  on  the  other  hand,  are  based  on  Taylor-series  expansions  of  y 
and  q about  y'  and  q',  truncated  after  the  first-order  term;  good  convergence 
depends  on  6q  and  by  being  small  relative  to  q and  y.  The  quadratic  term  of 
the  series  expansion  for  q was  approximately  taken  into  account  by  the 
following  procedure.  It  is  noted  that,  approximately, 

D,,,  = Pj,,  -H  W 

and  that 

6x 

= = , , (q, 

r'  2 .,2 

Yj 

/ 5 Cl  \2 

= Di  1 + ) + 

Qj 

while 

dP 

+ 1 -I-  = P + \Ni  + ^6yL 


(6.88) 

(6.89) 

(6.90) 

dP  aw  - 

-(-  ayp  -f-  — ^6y  (6.91) 

aYR  dy 


j + i 


, aD  ^ 
2 -f-  6y 


a y 


aD  ~ 
^ay 


'J 


Thus, 


1 -h 


5q 


r- 


Pj  Wj 


ayL 


SYl  + 


3Yr 


by, 


D, 


(6.92) 


On  the  assumption  that,  in  the  trailing  region,  R,^^changes  fromjDne  iteration 
to  the  next  primarily  as  the  result  of  changes  in  q,  rather  than  y (as  was 
often  found),  one  can  write 


/ P,  + w“ 

V 


1 


(6.93) 


the  absolute  value  signs  being  put  in  as  a precaution  against  accidental 
negative  values  for  P -t-  W— recall  that  all  these  values  are  numerically  very 
small  in  the  trailing  region  of  flow. 
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6.7  Solution  of  the 
Linear  System  by  the 
Double-Sweep  Technique 


Thus,  on  the  grounds  that  the  quantity  dRJdq^  represents  the  coefficient  of 
SQl  in  the  matrix  of  coefficients,  it  is  modified  by  the  factor  f to  the  result. 


in  which 


Worthy  of  note,  f differs  but  little  from  unity.  Similarly, 


9Rm 

3qR 


(6.94) 


(6.95) 


(6.96) 


There  remains  to  set 

^ = -Vr  + <^yp(bL  - bp)  -H  D ^ -)  (6.97) 

ayp  ' Ch  ay  y ' 


In  the  advance  phase,  in  the  tip  cell. 


aR,^  _ W - D 

axp  axLR 


(6.98) 


At  each  iteration  the  linear  system  represented  by  equations  6.77  must  be 
solved.  The  matrix  of  coefficients  is  sparse,  in  that  non-zero  values  are 
clustered  around  the  diagonal.  The  following  efficient  method  of  solution  that 
utilizes  this  sparseness  is  based  upon  a scheme  of  Preissman,  presented  by 
Liggett  and  Cunge  (4). 


The  given  system  can  be  represented  schematically  by  the  pair  of  equations 


a -h  b rjk  + c + d 77k+i  -1-9  = 0 
P + q -1-  r ^k+1  + S r7k+i  -I-  t = 0 


(6.99) 


one  such  pair  for  each  cell  k = 1,2,  . . .,  N.  Here  the  and  Tjk  represent  the 
corrections  6yk,  Spk  to  be  made  at  station  k in  the  given  iteration;  the  a,  b,  c,  d 
represent  the  appropriate  partial  derivatives  of  Rc,  while  g is  R^  itself.  The  p, 
q,  r,  s represent  the  partial  derivatives  of  R^,  and  t is  R^. 


Auxiliary  coefficients  Ek  and  Fk  are  defined  such  that 

»7k  = Ek  ^k  + Fk 


(6.100) 
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whatever  the  values  I,,,  prove  to  be.  Of  course  equation  6.100  does  not 
define  Ek  and  Fk  uniquely,  since  an  infinite  variety  of  value-pairs  will  satisfy  it 
for  any  given  pair  of  solution  values  rjk,  $k-  However,  any  pair  of  values  Ek,  Fk 
which  satisfy  equation  6.100  will  do.  For  example,  for  k = 1,  a correct  pair  is 


(6.101) 


in  which  is  known  from  the  upstream  boundary  condition. 

A scheme  is  sought  whereby  Ek+i  and  Fk+i  can  be  determined  in  terms  of  Ek 
and  Fk.  For  then,  all  E and  F could  be  determined  sequentially  for  k = 2,  3, 

. . .,  N -H  1 in  the  so-called  forward  sweep.  With  known  from  the  down- 
stream boundary  condition,  equation  6.100  would  yield  tjn  + i.  Furthermore, 
substitution  of  Tjk  from  equation  6.100  into  either  of  equations  6.99  would 
yield  an  expression  for  in  terms  of  rjk+i-  This  then  would  be  used  to 
obtain  from  + Equation  6.100  would  then  yield  tjn,  and  so  on,  with 

k = N-1,N-2,  ...,1,in  execution  of  the  backward  sweep. 

The  desired  relations  for  Ek+i  and  Fk+i  in  terms  of  Ek  and  Fk  are  readily,  if 
tediously,  found  by  substituting  rjk  from  equation  6.100  into  each  of  equations 
6.99,  then  solving  each  equation  for  and  equating  the  results.  Collection  of 
terms  yields  an  equation  of  the  form 


(6.102) 


in  which 


r Gi  - c Gg 


d G3  — s Gi 

G4  Gi  - Gg  G3 

d G3  - s Gi 


(6.103) 


with 


Gi  = a -f  b E| 


(6.104) 


G2  = g + b F| 


(6.105) 


G3  = p -I-  q E| 


(6.106) 


G4  = t -H  q F, 


(6.107) 
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In  the  backward  sweep 


= U + Wn,,,  + w 

Vk  = Ek  + Fk 


with 


(6.108) 


(6.109) 

(6.110) 

(6.111) 


6.8  Nondimensional 
Forms  of  the  Governing 
Equations 


The  program  works  in  a dimensioniess  mode.  This  aiiows  convenient  entry  of 
nondimensional  data  if  a series  of  generaiized  runs  are  being  made.  When 
dimensioned  input  data  are  used,  it  is  stiii  an  advantage  to  work  in  dimen- 
sioniess mode,  for  the  dimensionless  variables  are  also  normalized  in  some 
sense  as  weli;  that  is,  depths  and  discharges  in  the  stream  generaiiy  run  in 
the  range  0 to  1;  dimensionless  bottom  slope  in  sioping  borders  is  aiways 
unity,  and  so  forth. 


Dimensioniess  variabies  are  ratios  of  the  dimensioned  variabies  to  some  per- 
tinent characteristic  vaiue  having  the  same  dimensions.  Dimensioniess  varia- 
bles (starred)  are  defined  as  foilows 


Here  Yq  is  a characteristic  depth,  Qq  is  a characteristic  discharge  (per  unit 
width),  Xo  is  a characteristic  iength,  and  Tq  is  a characteristic  time. 

In  sloping  borders,  the  characteristic  depth  is  taken  as  normal  depth  for  the 
given  infiow  discharge  and  average  bottom  siope  and  roughness 
(DMLMOD  = 1), 


Yo  = Vn  (din,  So,  n) 


(6.113) 
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(with,  say,  the  Manning  formula  in  force).  The  average  bottom  slope  Sg  is 
given  by  the  total  change  in  bottom  elevation  divided  by  border  length.  The 
characteristic  discharge  is  simply  the  initial  inflow  rate 


Qo  — din 

The  characteristic  distance  is  given  by 


and  the  characteristic  time,  by 


(6.114) 


(6.115) 


(6.116) 


In  horizontal  borders,  normal  depth  is  infinitely  great,  so  alternative  defini- 
tions are  necessary.  With  DMLMOD  = 2, 

Qo  = din  (6.117) 


as  before,  and 


and 


To  — tc 


Xo  Yo  = Qo  T 


0 


(6.118) 


(6.119) 


as  before,  a convenient  choice  for  the  characteristic  depth  is 


Yo 


Q3/5  T V5 
0 ' 0 


(6.120) 


in  which  the  Chezy  Cg  is  evaluated  at  Yo. 

With  either  set  of  choices  for  the  characteristic  variables,  the  formulation  of 
dimensionless  drag  turns  out  particularly  simple. 

Substitution  of  the  products  y = y*  Yo,  x = x*  Xg,  d = d*  Qo.  etc,  into  the 
equation  of  mass  conservation,  6.12,  results  in  an  equation  of  exactly  the  same 
form,  but  with  all  variables  but  the  weighting  factors  starred,  that  is,  dimen- 
sionless. Only  the  computation  of  z is  affected  by  the  transition  to  dimen- 
sionless variables.  For  with 


z = T + c (6.121) 


■'“Note  that  b in  this  context  is  not  the  same  as  the  b used 
elsewhere  to  represent  bottom  elevation. 
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it  follows  that 


z*  = K*  r*®  + B*  T*  + C* 


(6.122) 


in  which 


K* 


(6.123) 


B* 


(6.124) 


(6.125) 


K*,  B*,  and  C*  are  evidently  dinnensionless  infiltration  coefficients  character- 
istic of  the  given  irrigation.  Worthy  of  note,  the  numerical  values  of  these 
parameters  are  different  in  the  two  schemes  of  converting  variables  to 
nondimensional  form. 

Similar  substitution  into  the  equilibrium  equation,  6.42,  yields  expressions  for 
P and  W in  the  starred  variables  that  are  identical  in  form  to  those  with 
dimensioned  variables,  equations  6.24  and  6.25.  Only  in  the  case  of  sloping 
borders,  the  dimensionless  bottom-elevation  difference  is  bL*  - bp*  = Xr*  - Xl*; 
in  horizontal  borders,  of  course,  bL*  - bp*  = 0.  Equation  6.39  for,  now,  the 
dimensionless  drag  intensity  retains  its  form,  except  that  the  Chezy  Ch  is  now 
a dimensionless  version,  different  from  equation  6.31.  For  both  choices  of 
characteristic  variables,  the  dimensionless  Chezy/Manning  Ch  is  simply 


Ch* 


(6.126) 


as  can  be  readily  proved  by  introducing  both  sets  of  dimensionless  variables 
into  the  equilibrium  equation. 

The  dimensionless  inflow  rate,  for  constant  inflow  to  cut  off,  is  simply  unity. 

In  the  second  system,  t^o*  = 1 by  definition;  in  the  first  system,  it  is  a 
variable.  Field  end  is  marked  by  x*  = L*;  of  note,  the  numerical  value  of  L*  is 
different  in  the  two  systems,  for  a given  real,  dimensioned  L In  general,  the 
variables  of  interest  have  differing  numerical  values  for  the  same  irrigation, 
viewed  in  the  two  systems  of  nondimensionalization. 

In  the  first  system  (DMLMOD  = 1),  the  dimensionless  solution— controlling 
variables  (as  described  in  (8))  are,  with  uniform  bottom  slope  and  Manning 
resistance,  the  field  length  L*,  the  infiltration  parameters  a,  K*,  B*  C*,  and 
cut-off  time  too*.  In  the  second  system,  these  are  slope  So*  and  L*,  and  a,  K* 
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6.9  Differential  Forms 
of  the  Governing 
Equations 


For  comparison  with  other  studies  in  the  literature,  it  is  worthwhile  to  express 
the  integrated  equations  (integrated  over  cell  length  and  a time  increment)  of 
mass  conservation  and  equilibrium  in  differential  form,  as  partial  differential 
equations  in  x and  t. 


The  equation  of  mass  (volume)  conservation,  equation  6.12,  is  first  written  for 
a very  short  cell— short  enough  for  the  <j)  to  be  one-half— and  for  a short  time 
increment.  Then  it  is  divided  by  6Xlr  and  by  6t.  Then  6Xlr  and  6t  are  both  made 
to  approach  zero.  The  weighted  sums  then  assume  much  simpler  forms.  For 
example  0qL  -t-  (1  - 0)  qj  — qL,  etc.  Furthermore,  in  the  limit,  as  N — cx. 


6Xlr 


(Yl  + Zl) 


6Xlr 


(Vr  + Zr) 


1 dxA  d{y  + z) 



N dt  dx 


(6.127) 


approaches  zero  faster  than  (q,^  + qR)/6xLR,  and  the  “scooping  in’’  of  volume 
(the  second  term  in  the  right  of  each  of  equations  6.13  and  6.14),  stemming 
from  the  obliqueness  of  the  computational  grid,  becomes  negligibie.  The  result 
in  the  limit  is  the,  so  called,  equation  of  continuity, 

^ ^ ^ 0 (6.128) 

ax  at  at 


since 

Yl  + 0yR  Yr  = y ^ y.  etc.  (6.129) 

The  equilibrium  equation,  6.42,  written  for  a thin  slice,  following  division  by 
the  volume  of  the  slice  y 6xlr  and  passage  to  the  limit  as  6Xlr  — 0,  yields  the 
result 


A = So  - S,  (6.130) 

ax 

in  which  Sf  is  given  by  equation  6.28,  (eq.  6.129  was  used  in  obtaining 
eq.  6.130). 


/ 
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Part  II:  Program  Organization 


6.10  The  First  Time  Step 


6.11  Preparation  of  Trial 
Values  for  Each 
Successive  Time  Line 


Time  and  advance  are  set  to  zero  for  the  time  index  i = 1.  The  computation 
of  the  first  time  step  seeks  the  advance  and  profile  at  i = 2.  The  time  t(2)  = St^,  in 
which  5ti  is  the  size  of  the  first  time  step. 

Dimensionless  values  of  R.  and  R„  are  each  set  to  0.000001.  The  weight- 

'"max  '"max  ^ 

ing  factor  in  time,  6,  is  set  to  0.6. 

The  advance  in  the  first  time  step  is  computed,  with  the  entire  stream  length 
treated  as  a single  (tip)  cell,  and  appropriate  weighting  factors  used  in  the 
equations.  It  is  assumed  that  discharge  increases  according  to  a three- 
seventh  power  law  from  0 to  the  value  Pin  in  the  course  of  the  first  time  step; 
consequently  the  volume  of  inflow  in  the  first  time  step  is  0.7  q^  5ti.  This  first 
step  of  calculation  is  always  handled  in  fully  nonlinear  fashion,  regardless  of 
whether  subsequent-time-step  calculations  are  locally  linearized  or  not.  Once 
the  stream  length  and  upstream  depth  have  been  determined  for  the  first  time 
step,  the  length  is  subdivided  into  N cells. The  depths  y^  at  the  juncture 
points  k = 2,  3,  4,  . . .,  N,  are  found  by  means  of  equation  6.52  with  /3  taken 
from  equation  6.51  (the  total  length  is  substituted  for  [Xg  - x^]).  Infiltration 
depths  Zk  are  found  from  equation  6.61  (with  the  same  substitution  of  total 
stream  length  for  8x^).  The  discharges  Pk  are  then  given  by  the  formula 

qk  = qLB*^^^^ ^ (6.131) 

Ylb  + Zlb 


The  program  is  then  sent  to  section  6.18,  End  of  Time  Step. 

In  equation  6.131,' the  subscript  LB  refers,  as  before,  to  conditions  at  the  left 
(upstream)  stream  boundary;  RB  refers  to  the  downstream  (right)  stream  boun- 
dary. It  follows  that  k = 1 corresponds  to  LB,  while  k = N -i-  1 corresponds 
to  RB.  During  advance,  Xg  refers  to  Xrb-  The  subscripts  L and  R refer  to  k and 
k -I-  1,  respectively. 

For  successive  time  steps,  the  procedure  is  to  increase  the  time  index  i by 
one,  and  increment  the  time  by  6t.  The  time  increment  varies  in  accordance 
with  one  or  another  scheme,  defined  by  the  choice  of  DTMOD  (sec.  3.9).  In 
addition,  it  is  convenient  to  compute  the  ultimate  total  volume  of  inflow  as 
Vq  = Pin  tco,  in  which  q^  is  a constant  inflow  rate,  and  tco  is  the  cut-off  time. 
On  the  other  hand,  the  volume  of  inflow  during  the  first  time  step  is  set  to 


^■'See  section  3.9  (DTMOD)  for  possible  variations  on  this 
scheme. 
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0-7  • Clin  • (sec.  6.10).  Further,  the  cut  off  is  smeared  over  one  time  step,  so 
that  the  increment  of  inflow  volume  pertaining  to  the  step  that  contains  t^o  is 
(1  - 0)qin6t.  To  reconcile  all  these  factors,  one  time  step  is  made  just  the 
right  size  to  ensure  that  a time  line  falls  precisely  on  the  precutoff  value 
given  by 


tpco  = U - (1  - 0)6t3,d  + 0.3  6t,  (6.132) 

Then,  t^o  vvill  fall  at  just  that  fraction  of  the  standard  time  step  after  tp^o  that 
will  allow  the  inflow  volume  computed  by  increments  over  all  the  time  steps 
to  equal  . tpo- 

In  any  event,  if  t|  exceeds  the  specified  t,^ax,  fhe  computations  are  ended  in 
one  of  the  two  normal  modes  of  termination.  The  other,  standard,  normal  ter- 
mination follows  post-irrigation  analysis,  when  the  computed  irrigation  has 
been  completed. 

With  computations  continuing,  first  guesses  are  made  for  Xlb,  Xrb,  and  x^,  yk, 

Pk  for  all  k.  The  left  and  right  boundaries  are  assumed  the  same  as  on  the 
previous  time  line,  except  during  advance;  then,  for  a first  guess,  the  right 
boundary  is  assumed  to  advance  at  the  same  speed  as  in  the  previous  time 
step.  The  Xk  are  found  by  subdividing  the  total  computational  stream  length 
Xrb  - Xlb  into  N parts.  The  yk  and  Pk  are  assumed  equal  to  their  values  yknownk- 
Pknowrik  on  the  previous  time  line.  Bottom  elevation  and  depths  of  infiltration 
are  computed  for  all  k.  The  value  of  discharge  at  the  upstream  boundary  is 
computed  and  the  correction  6qi  = r?i  = q^  - Pknown^  to  the  first  guess  for  dis- 
charge at  LB  is  obtained. 


6.12  Overview  of  the 
Iteration  Cycle  for 
Correcting  Trial  Values 


At  this  point  the  iteration  cycle  is  entered.  If  the  locally  linearized  mode  of 
solution  has  been  chosen,  the  cycle  is  traversed  just  once.  If  the  nonlinear, 
iterative  mode  is  in  force,  the  cycle  is  traversed  repeatedly,  the  iteration  index 
j increased  by  1 each  time,  until  Rp  and  in  every  cell  are  reduced  below 
Rp,^^^  and  Rm^ax’  respectively,  or  j exceeds  jmax-  In  the  last  case,  the  time  step 
is  cut  in  half  and  the  calculations  repeated,  and  so  on,  for  a maximum  of 
eight  reductions.  If  convergence  is  still  not  achieved,  the  error  message, 

“NO  CONVERGENCE  IN  CELLS”  is  printed,  and  an  abnormal  exit  (controlied^^ 
abort)  is  made. 


. .as  opposed  to  uncontrolled  aborts,  defined  here,  as 
aborts  initiated  by  computer-system  software.  Any 
diagnostic  information  printed  is  then  controlled  solely  by 
that  software. 
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6.13  Cell  Computations 
and  Solution  of  the 
Resulting  Set  of  Linear 
Equations 


If  Rc  or  Rm  in  any  cell  (eqs.  6.68,  6.69)  exceed  the  permitted  maximums,  note 
is  taken  of  that  fact.^^  The  partial  derivatives  (eqs.  6.79-98)  comprising  the 
coefficients  of  the  corrections  r?,,,  in  equations  6.77  are  computed  and  the 
equations  solved  by  the  double-sweep  method  outlined  in  section  6.7.  The 
corrections  for  k = 1,  2,  . . .,  N -i-  1 are  added  to  the  current  (j""")  approxima' 
tions  of  the  variables. 


6.14  Elimination  of 
Computed  Negative  Cell 
Lengths:  Restructuring 
System  of  Nodes 


If  the  stream  is  in  the  advance  phase,  every  new  profile  is  checked  immedi- 
ately for  negative  cell  lengths.  These  can  arise  because  the  node  distribution 
is  based  on  the  first  guess  for  the  advance  increment.  With  closely  spaced 
nodes  near  the  leading  edge,  a negative  correction  to  the  advance  increment 
can  place  the  new  right  boundary  upstream  from  an  interior  node  point.’'* 


The  locations  of  the  interior  node  points  are  not  generally  changed  in 
response  to  the  iterative  adjustment  of  the  stream  boundaries,  because  this 
would  add  the  location  of  the  right  stream  boundary  as  unknown  in  the  equa- 
tions for  each  cell,  instead  of  for  just  the  downstream-boundary  cell.  This,  in 
turn,  would  change  the  appearance  of  the  matrix  of  coefficients  for  the 
system  of  linear  correction  equations.  Instead  of  restricting  the  non-zero  coef- 
ficients to  a narrow  band  about  the  diagonal,  the  system  would  also  have  an 
additional  vertical  stripe  of  coefficients  applying  to  the  unknown  stream 
length.  This  would  require  a major  change  in  the  method  of  solution  of  the 
equation  set. 


Instead,  the  internal  nodes  are  left  fixed  in  the  hope  that  only  the  length  of 
the  last  node  will  change  in  consequence  of  corrections  to  the  advance  incre- 
ment. Only  in  the  event  that  a negative  correction  exceeds  the  length  of  the 
right  boundary  cell  does  it  become  necessary  to  restructure  the  system  of 
nodes.  Then,  a new  node  set  is  constructed,  in  accord  with  the  given  value  of 
RDX,’^  but  based  on  the  corrected  stream  length.  When  this  occurs,  a recom- 
putation for  the  time  step  is  made,  that  is,  the  iterations  are  restarted  from 
the  beginning,  with  the  original  trial  values  of  depth,  discharge,  and  advance 
increment.  Only  the  iteration  counter  continues  to  increment,  so  that  the  total 
number  of  iterations  is  always  noted  and  does  not  exceed  the  maximum,  j,^ax- 
Thus,  even  with  the  locally  linearized-solution  mode  (LINMOD  = 1),  the 
number  of  iterations  can,  on  occasion,  exceed  one. 


''^RFLAG  is  set  to  one,  in  contrast  to  the  zero  value  set  at 
the  start  of  the  iterations. 

. . .either  because  the  ratio  RDX  is  small,  the  number  of 
cells  NSTD  is  large,  or  the  time  step  DT  is  large  (see  Input, 
line  10,  sec.  3.10). 

■■^Input,  line  10,  sec.  3.10. 
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Once  the  stream  length  is  seen  to  consist  of  positive  cell  lengths  only,  all  of 
the  new  depths  are  scanned  for  negative  values,  for  the  program  operates 
under  the  general  rule  that  no  surface-stream  computations  will  be  performed 
with  negative  depths.  Apart  from  the  physical  impossibility  of  a negative 
depth,  certain  terms  in  the  governing  equations,  notably  the  drag  term,  have 
no  mathematical  meaning  either,  with  a negative  or  zero  depth,  and  any 
attempt  to  calculate  drag  with  a non-positive  depth  will  result  in  an  unsched- 
uled, system-controlled^®  abort.  Computed  depths,  then,  less  than  zero  simply 
are  not  allowed  in  the  surface  stream. 

Apparent  negative  values  can  arise  from  either  of  two  causes,  one  completely 
normal,  the  other  evidence  of  a less-than-ideal  computation.  As  regards  the 
first  cause,  once  inflow  is  cut  off,  depths  at  the  rear  of  the  surface  stream 
continually  decrease  because  of  infiltration  and,  usually  also,  because  of 
bulk  flow  in  the  downstream  direction.  In  a stream  with  incipient  front-end 
recession,  the  depths  near  the  right  stream  boundary  also  decrease.  Since 
each  computed  approximation  to  a variable  value  is  based  on  equations 
linearized  about  the  current  value  of  the  variable,  and  these  are  applied  over  a 
full  time  step  of  some  given  size,  if  the  true  result  calls  for  a zero  depth  at 
some  time  within  the  duration  of  the  time  step,  the  computed  depth  at  the 
end  of  the  time  step  will  be  negative.  This  is  a normal  computation  during 
recession,  and  the  time  of  passage  of  the  recession  edge  past  a given 
x-station  is  estimated  by  linear  interpolation  between  the  computed  depths  at 
that  station  at  the  beginning  and  end  of  the  time  step.  In  a completely  normal 
calculation,  if  several  stations  experience  computed  depths  passing  to 
negative  values  in  a given  time  step,  the  station  with  the  largest^®  computed 
negative  depth  lies  at  the  boundary  of  the  stream,  and  the  neighboring  sta- 
tions exhibit  computed  negative  depths  progressively  smaller^®  towards  the 
interior  of  the  stream  (or,  of  course,  positive  depths  progressively  larger). 

If  the  program  is  working  well,  in  the  sense  that  all  functions  of  interest  are 
changing  smoothly  and  gradually  with  respect  to  time  and  distance,  no  nega- 
tive depths  will  be  computed  in  the  interior  of  the  profile  unless  all  depths  to 
the  exterior  of  that  point  are  also  negative,  as  described  above.  An  imperfect 
calculation,  on  the  other  hand,  can  lead  to  the  development  of  saw-tooth  pro- 
files, in  which  some  depth  peaks  are  positive  and  some  troughs  negative.  In 
particular,  the  drag  term  in  the  force-equilibrium  equation  is  very  sensitive  to 
small  errors  in  discharge  or  depth  when  the  depth  is  very  low.  These  small 
errors  can  lead  to  large  errors  in  drag,  with  subsequent  computation  of  saw- 


6.15  Elimination  of 
Computed  Negative 
Stream  Depth: 
Recession 


. .as  opposed  to  program-controlled. 

■'^Zero  depths  are  allowed  only  at  the  stream  boundaries,  for 
drag  is  not  computed  there. 

. . in  the  absolute  sense.  . . 
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tooth  profiles.  Furthermore,  use  of  imprecise  weighting  factors  when  cell 
lengths  are  large  can  also  lead  to  saw-tooth  profiles.  If  negative  depths  arise 
when  all  the  surface  depths  are  quite  low  near  the  end  of  the  irrigation, 
reduction  of  the  computational  stream  profile  to  those  regions  that  are  com- 
prised of  positive  depths  only,  is  a well  justified  means  of  getting  useful  final 
results  from  a troublesome  calculation.  If  the  saw  teeth  are  of  high 
amplitude,  occurring  in  a region  of  large  depth,  however,  then  significant 
volumes  of  water  will  be  lost  in  the  calculation,  and  recession,  either  from 
the  rear  or  from  the  front,  will  be  computed  as  occurring  too  rapidly.  This  con- 
dition is  signaled  by  a large,  positive  relative  volume  error  (see  post-irrigation 
analysis,  sec.  5.5),  well  over  0.01;  this  circumstance  alerts  the  user  to  the 
problem.  The  procedure  for  eliminating  those  portions  of  the  profile  that  con- 
tain negative  depths  is  outlined  next. 

The  entire  profile  between  Xlb  and  Xrb  is  scanned  for  negative  depths  and  the 
location  of  each  occurrence  is  noted.  If  none  is  found,  the  program  proceeds 
to  6.16.  If  any  are  found,  the  largest  length  of  stream  either  between  a boun- 
dary and  the  nearest  location  of  negative  depth  or  between  two  locations  of 
negative  depths  is  selected  as  a basis  for  the  new  computational  stream 
length.  These  new  boundaries  are  considered  tentative,  pending  examination 
of  the  depths  there. 

If  either  boundary  has  been  stepped  in  toward  the  interior,  the  entire  time 
step  will  be  recomputed  on  the  basis  of  the  new  stream  length.  This  length 
now  applies  to  both  the  accepted  profiles  at  the  beginning  of  the  time  step 
and  the  unknown  profile  at  the  end  of  the  time  step.^®  In  other  words,  only  a 
portion  of  the  “old”  stream— albeit  the  most  significant  portion— is  used  to 
advance  the  calculations  to  the  new  time  line.  On  the  other  hand,  if  the 
largest  stream  length  following  the  step-in  process  is  less  than  1 percent  of 
field  length,  the  irrigation  is  presumed  over.  The  trailing  recession  edge  is 
assumed  to  cross  the  new  upstream  boundary  at  the  current  time.  Front-end 
recession  is  assumed  to  cross  the  new  downstream  boundary  also  at  the  cur- 
rent time.  If  such  stepping  in  of  a computational  boundary  constitutes  the 
first  indication  of  rear-end  or  front-end  recession,  recession  from  the 
corresponding  old  stream  boundary  is  assumed  to  occur  arbitrarily  a half 
time  step  earlier.  Once  recession  stations  and  times  are  noted,  program  con- 
trol passes  to  post-irrigation  analysis  (sec.  5.5) 

With  a stream  of  working  length,  the  depth  at  the  tentative  upstream  boun- 
dary is  checked  first.  If  positive,  the  location  is  accepted  as  the  new  boundary 
and  the  tentative  downstream  boundary  is  viewed  next.  But  if  the  depth  there 


''^The  iteration  counter  continues  to  be  incremented,  rather 
than  set  back  to  one. 


is  negative,  recession  has  occurred  from  the  old  upstream  boundary.  If  the 
depth  at  the  old  boundary  location^®  has  turned  negative,  the  time  of  reces- 
sion there  is  determined  by  linear  interpolation  between  the  depths  computed 
there  at  the  two  times.  If  the  depth  is  still  positive,  signifying  a saw-tooth  pro- 
file, recession  from  the  old  boundary  is  still  assumed  to  have  occurred, 
arbitrarily  one-half  of  a time  step  earlier  than  the  current  time. 

Then,  a location  with  a small  positive  depth  in  the  new  profile  is  sought  for 
the  new  computational  boundary.  The  profile  downstream  from  the  tentative 
boundary  location  is  scanned  for  a depth  equal  to  or  greater  than  a preset 
small  value,  yrec-^^  If  all  depths  are  less  than  y,ec,  the  irrigation  is  assumed 
over.  Any  remaining  surface  water  is  added  to  the  infiltrated  depths  by  com- 
puting recession  times  commensurate  with  the  additional  time  necessary  to 
infiltrate  those  small  depths.  The  program  then  passes  to  post-irrigation 
analysis. 

But  if  a depth  greater  than  y,ec  is  found  at  an  interior  node,  the  new  computa- 
tional boundary  is  placed  at  the  location  nearest  the  tentative  boundary  at 
which  a depth  equal  to  y^ec  exists,  as  determined  by  linear  interpolation 
between  nodes.  In  preparation  for  a recomputation  of  the  same  time  step 
with  the  new  boundary,  the  upstream  portion  of  the  computational  stream  as 
it  existed  at  the  beginning  of  the  time  step  is  also  truncated  to  this  same  new 
x-boundary.  Only  the  portions  interior  to  this  point  will  enter  into  the  volume- 
balance  equations.  The  depth  and  discharge  on  this  previous  time  line  are 
found  at  this  x by  linear  interpolation  between  nodes;  the  depth  of  infiltration 
there  is  found  by  applying  the  infiltration  function  to  the  appropriate  infiltra- 
tion time.  Cell  volumes  for  the  new  upstream-boundary  cell  are  found  by 
using  the  appropriate  weighted  averages  of  nodal  surface  and  subsurface 
depth  values. 

In  the  stepping-in  process,  the  total  number  of  cells  in  the  stream  may 
change.  If  the  boundary  is  moved  more  than  one  cell  length  from  its  previous 
position  on  the  old  time  line,  at  least  one  cell  will  be  lost  from  consideration. 
The  node  number  at  the  left  boundary  is  taken  to  be  the  number  of  the 
nearest  node,  on  the  old  time  line,  upstream  from  the  point  where  y = y,eo  on 
the  new  time  line.  This  insures  that  the  number  of  cells  on  the  new  and  old 
time  lines  are  the  same.  (After  the  downstream  boundary  has  been  dealt  with, 
and  prior  to  recomputation  of  the  time  step,  the  nodes  are  renumbered,  so 
that  the  first  node  is  always  indexed  one.) 


^°The  old  and  new  boundaries  can  be  the  same. 

^''(dimensionless)  y,ec  = 0.01. 
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The  same  considerations  are  now  applied  to  the  tentative  downstream  boun- 
dary. If  the  depth  there  is  nonnegative,  the  instructions  described  in  the  rest 
of  this  paragraph  are  bypassed.  In  the  event  of  a negative  depth,  the  boun- 
dary is  stepped  in,  following  calculation  of  a recession  time  for  the  old 
downstream  boundary  location,  to  the  nearest  location  where  a surface  depth 
equals  y^ec-  If  none  is  found,  the  irrigation  is  over;  the  remaining  surface 
depths  are  added  to  the  infiltration  profile,  and  the  program  passes  to  post- 
irrigation analysis  (sec.  5.5).  If  a new  boundary  is  located,  however,  the  old 
computational  stream  length  is  again  truncated.  If  the  resultant  stream  length 
is  less  than  1 percent  of  the  field  length,  the  recession  edges  are  judged  to 
pass  the  new  boundaries  at  the  current  time,  and  the  irrigation  is  deemed 
over;  post-irrigation  analysis  commences.  Otherwise,  the  dimensions  and 
volumes  of  the  downstream-boundary  cells,  on  the  previous  and  on  the  cur- 
rent time  lines,  are  recalculated  as  for  the  upstream-most  cell,  prior  to  the 
recomputation  of  the  entire  time  step. 


In  this  way,  even  if  a considerable  shift  in  boundary  location  occurs  in  the 
course  of  a single  time  step,  the  computational  grid  never  becomes  very 
oblique.  A cell  boundary  that  is  moved  a great  distance  in  the  course  of  a 
time  step  can  lead  to  considerable  error  in  the  cell  mass-balance  computa- 
tions, if  the  infiltration  profile  is  relatively  deep  and  nonuniform. 

Since  whole  cells  may  have  been  lopped  off  the  stream  in  the  step-in 
process,  a renumbering  of  the  nodes  now  takes  place.  The  new  total  number 
of  cells  is  noted.  The  new  upstream-boundary  node  is  numbered  one,  and  the 
remaining  nodes  are  numbered  sequentially  progressing  downstream.  The 
arrays  describing  conditions  at  each  node,  new  and  old  x,  bottom  elevation, 
surface  depth,  discharge,  and  infiltration  time  and  depth,  are  shifted  in 
accordance  with  the  new  numbering. 


At  this  point,  that  portion  of  the  infiltrated  volume  that  lies  at  the  begin- 
ning of  the  time  step  outside  the  current  boundaries  of  the  stream  is 
calculated  and  stored  for  subsequent  use  by  the  volume-balance  check  that 
is  performed  for  diagnostic  purposes  at  every  time  step. 


Next,  the  new  conditions  at  the  right  boundary  are  investigated. 


6.16  Investigation  of 
Potential  Changes  in 
Downstream  Boundary 
Conditions 


Apart  from  the  possibility  of  front-end  recession  signaled  by  the  computation 
of  negative  depths,  which  has  been  treated  in  the  previous  section,  a number 
of  other  possibilities  arise  at  the  downstream  boundary  of  the  surface  stream, 
and  the  program  must  branch  accordingly.  The  first  branching  is  determined 
by  conditions  at  the  beginning  of  the  time  step.  Each  of  these  paths  then 
branches  further,  depending  upon  what  has  developed  during  the  time  step. 
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6.16.1  Advance  at  Beginning  of  Time  Step 

A stream  that  was  advancing  at  the  beginning  of  the  time  step  will  be  com- 
puted for  the  end  of  the  time  step  as  (1)  continuing  advance,  (2)  exhibiting 
negative  advance,  that  is,  spontaneously  receding  from  the  front  end,  (sec. 
6.16.1.1)  or  (3)  crossing  the  field  boundary  (sec.  6.16.1.2). 

When  a positive  advance  increment  is  computed,  the  stream  is  ready  for  the 
next  iteration,  or  the  next  time  step,^^  or  a change  in  boundary  condition 
because  field  end  has  been  crossed. 

6.16.1.1  Front-End  Recession  at  End  of  Time  Step 

A computed  negative  advance  increment  requires  a recomputation  of  the  time 
step.  If  the  nonlinear  approach  has  been  chosen,  the  spontaneous  front-end 
recession  is  not  treated  as  such,  that  is,  recomputation  does  not  commence 
until  the  errors  in  mass  and  force  balances  have  been  reduced  below  the 
required  tolerances  in  the  iterations,  thus  confirming  the  computation  of 
negative  advance.  With  the  linear  approach,  the  first  indication  is  sufficient. 

In  any  event,  the  recomputation  is  performed  on  the  basis  of  front-end  reces- 
sion from  a boundary  located  at  the  point  of  maximum  advance  (that  is,  the 
location  of  the  right  boundary  at  the  beginning  of  the  time  step).^^ 

6.16.1.2  Stream  Encounters  Field  End:  Runoff  or  Ponding  at  End  of  Time  Step 

When  an  advance  computation  shows  the  stream  front  crossing  field  end,  in 
fact,  advance  was  completed  at  some  point  during  the  current  time  step,  and 
either  runoff  or  ponding  is  occurring  for  the  rest  of  the  step.  This  calls  for  a 
recomputation  of  the  time  step  with  the  new  boundary  conditions  applied  to 
the  current  time  line.  In  the  nonlinear  approach,  the  iterations  are  allowed  to 
continue  normally  to  confirm  stream  arrival  “past”  field  end.  The  time  of 
arrival  of  the  wave  at  field  end  is  found  by  assuming  a constant  rate  of 
advance  over  the  time  step  equal  to  that  computed  for  the  time  step  in  the 
“absence”  of  a downstream  boundary. 

The  imposition  of  a new  right-hand  boundary,  at  x = L,  calls  for  a reposition- 
ing of  the  nodes.^'’  These  are  placed  in  accord  with  the  number  of  cells 
remaining  and  with  the  given  ratio  r^x  of  successive  cell  lengths. 

22.  . .whichever  is  appropriate,  depending  upon  (1)  whether  a 
recomputation  has  been  ordered  because  of  boundary  step- 
in,  (2)  whether  the  soiution  technique  has  been  chosen  iocai- 
iy  iinearized  or  fuiiy  noniinear,  and  (3)  whether  the  errors  in 
mass  and  force  baiances  are  smaii  enough  (sec.  6.17). 

23|f  this  negative  advance  increment  arises  prior  to  infiow  cut 
off,  the  computation  is  judged  a faiiure,  a message  is 
printed  to  that  effect  and  a controiied  abort  is  executed. 

2'*Without  repositioning,  it  is  possibie,  for  exampie,  for  the 
new  boundary  iocation  to  iie  upstream  of  the  downstream- 
most  ceii(s). 
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For  runoff,  the  first  guesses  for  the  new  nodal  values  of  depth  and  discharge 
are  determined  by  interpolation  among  the  last-computed  values.  Only  the 
downstream-boundary  value  of  depth  is  set  to  zero.  The  stream  is  now  ready 
for  recomputation. 

For  the  ponding  case,  the  first  guess  is  arrived  at  by  assuming  that  the  sur- 
face and  subsurface  volumes  for  that  portion  of  the  stream  computed  as 
lying  beyond  field  end  (so-to-speak,  “spilled”  volume)  are  retained  upstream 
from  the  end  barrier  in  a standing  pool.  First  the  spilled  volume  Is  computed; 
then  the  nodes  are  repositioned  in  the  new  stream  length,  and  tentative  nodal 
values  of  depth  and  discharge  are  found  by  interpolation  among  the  last- 
computed  values.  The  dimensions  of  the  pool  are  found  by  assuming  its 
surface  level,  except  in  a transition  cell  in  which  the  water-surface  elevation 
varies  linearly  from  the  interpolated  value  on  the  upstream  side  to  pool  level 
on  its  downstream  side;  pool  volume,  of  course,  is  known.  All  first-guess 
discharges  in  the  pool  are  assumed  zero.  The  stream  is  now  ready  for 
recomputation. 

6.16.2  Runoff  at  Beginning  of  Time  Step 

If  the  right  boundary  was  subject  to  runoff  at  the  beginning  of  the  time  step, 
and  the  computed  discharge  there  at  the  end  of  the  time  step  is  positive,  this 
iteration  is  assumed  completed,  and  the  program  passes  to  the  instructions 
in  section  6.17. 

If,  however,  the  discharge  turns  negative  there,  a recomputation  is  indicated, 
with  new  boundary  conditions  reflecting  front-end  recession.  First,  though, 
with  a nonlinear  computation,  the  iterations  are  allowed  to  run  their  course, 
to  confirm  the  negative  runoff.  A further  test,  particularly  pertinent  to  linear 
calculations,  is  applied  to  help  ensure  that  the  negative  runoff  computed  truly 
signals  the  start  of  front-end  recession,  rather  than  an  indication  of  saw-tooth 
behavior  of  a small-valued  computed  runoff  function.  In  incipient  front-end 
recession,  the  stream  profile  is  concave  upwards,  while  the  profile  upstream 
from  a discharging  overfall  Is  concave  downward.  The  last  four  nodal  values 
of  depth  are  examined  to  determine  the  sense  of  concavity  of  the  profile.  If 
the  profile  appears  to  be  concave  upward,  it  is  taken  as  bona  fide  indication 
of  front-end  recession.  With  a concave-downward  profile,  the  negative  runoff 
is  judged  spurious,  and  the  runoff  phase  is  allowed  to  continue  in  the  expec- 
tation that  positive  saw-tooth  peaks  will  be  computed  shortly,  or  that  soon 
the  depth  profile  will  turn  concave  upward  and  allow  computed  front-end 
recession  to  commence. 


6.16.3  Ponding  or  Front-End  Recession  at  Beginning  of  Time  Step 


6.17  Branching  at  the 
End  of  the  Iteration  Cycle 


6.18  Calculations  at  the 
End  of  the  Time  Step 


For  a stream  undergoing  ponding,  or  in  front-end  recession  at  the  start  of  the 
time  step,  the  iteration  is  complete  at  this  point. 

At  this  point  the  program  notes  whether  the  recompute  flag  has  been  raised. 

If  so,  the  iteration  cycle  is  repeated  at  once.  Otherwise,  it  is  repeated  only  if 
the  nonlinear  mode  of  solution  is  in  effect,  and  if  Rg  or  have  exceeded  the 
allowable  maximum  in  any  cell. 

Once  the  computed  stream  geometry  and  discharge  distribution  for  the  cur- 
rent time  line  have  been  accepted,  certain  information  pertinent  to  the  profile 
is  stored.  The  depths  at  both  ends  of  the  profile  are  compared  with  the  max- 
imums  attained  there,  up  to  the  current  time  line,  and  form  new  maximums  if 
the  old  ones  were  exceeded.  Further  the  maximum  depth  in  the  profile  and  its 
location  are  noted  and  compared  with  the  global  maximum  for  the  irrigation, 
replacing  the  last  if  it  is  greater. 


A measure  ejs  of  profile  irregularity,  the  lack  of  smoothness  is  computed 


©ss  = (yk+2  + Yk  - 2yk  + i)  - (yk+1  + yk-1  - sy^)  (6.133) 

for  all  k between  2 and  N - 1;  the  maximum  in  the  profile,  TOOTH,  is  deter- 
mined and  compared  with  a global  maximum. 


Current  surface  and  subsurface  volumes,  Vy  and  respectively,  are  computed, 
and  their  sum  is  compared  to  the  difference  between  the  total  volumes  of 
inflow  Vq  and  runoff  V,o  in  a volume  balance.  The  volume  error 

Verr  = Vq  - Vy  - V,  - V,,  (6.134) 


is  divided  by  the  total  inflow  volume  to  obtain  the  relative  volume  error 


PC„ 


(6.135) 


Total  inflow  and  runoff  volumes  (for  the  entire  extant  period  of  irrigation)  are 
incremented. 


Profiles  and/or  an  abbreviated  current  description  of  the  stream  are  printed,  if 
called  for  by  the  extant  value  of  IDIAG  (sec.  3.11).  Profiles  are  plotted  at  this 
time  upon  demand  (sec.  3.12). 

The  values  Xknown(k),  bknown(k),  yknown(k),  qknown(k),  replaced 

by  x(k),  b(k),  y(k),  q(k),  6Vy(k),  6V,(k),  respectively  for  all  k = 1,2,  . . .,  N -t-  1. 

The  program  returns  to  section  6.11,  Preparation  of  Trial  Values  for  Each  Suc- 
cessive Time  Line. 
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7.  The  Kinematic-Wave  Model 


This  chapter  details  the  theory  and  program  logic  outlined  in  sections  5.3  and 
5.4. 


7.1  Introduction:  Key 
Assumption— Relatively 
Large  Bottom  Slope 


The  kinematic-wave  model  is  a numerical  solution  of  a family  of  ordinary  dif- 
ferential equations.  In  this  introductory  section,  the  governing  partial  differen- 
tial equation  leading  to  this  family  is  derived.  Reduction  to  the  ordinary  dif- 
ferential equations  and  their  solution  is  relegated  to  subsequent  sections. 


In  chapter  6,  the  unsteady  flow  in  the  surface  stream  was  analyzed  by  break- 
ing it  up  into  computational  cells,  or  slices.  The  depths  and  discharges  at  the 
cell  boundaries  and  the  total  stream  length  were  determined  in  a series  of 
time  steps  by  the  solution  of  algebraic  equations  expressing,  at  each  time 
step,  mass  conservation  and  equilibrium  of  forces  for  the  water  in  each  cell. 
The  algebraic  equations  were,  in  fact,  numerical  approximations  of  time  and 
space  integrals  entering  into  the  exact  expressions  of  mass  conservation  and 
force  equilibrium.  So,  in  general,  one  could  say  that  the  surface  flow  was 
analyzed  through  integral  relationships  applied  to  finite-sized  slices  of  the 
surface  stream  for  finite-sized  time  steps. 

Alternately,  the  problem  could  be  viewed  mathematically  as  expressed  in 
terms  of  partial  differential  equations.  These  are  readily  derived  from  equa- 
tions 6.12  and  6.23  (as  shown  in  sec.  6.9)  and  comprise  the  so-called  equa- 
tions of  continuity 


^ ^ ^ = 0 

dx  at  at 


(7.1) 


and  motion 


= S.  - S,  (7.2) 

Here,  x and  t are  distance  and  time  coordinates,  respectively,  q,  y,  and  z are 
local,  instantaneous  discharge,  surface  depth,  and  infiltrated  depth,  respec- 
tively, So  is  the  local  bottom  slope,  and  Sf  is  the  friction  slope  at  the  given  x 
and  t (defined  by  eq.  6.28.) 

The  forces  on  the  water  in  each  cell  were  assumed  to  arise  from  three 
sources;  namely,  the  component  of  the  weight  in  the  direction  of  flow,  the 
drag  of  vegetation  and  bottom  directed  opposite  to  the  flow,  and  the  dif- 
ference in  pressures  exerted  at  the  cell  boundaries  by  the  water  in  the 
neighboring  slices  upstream  and  downstream.  Now  for  those  flows  in  which 
the  weight  component  proves  much  larger  than  the  pressure  difference,  a 
major  mathematical  simplification  becomes  tenable.  For  then,  equation  7.2  is 
approximated  by  the  equality. 
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So  — S( 


(7.3) 


which  is,  in  fact,  an  algebraic  or  transcendental  relation  between  depth  and 
discharge  at  every  point  in  the  flow  (eq.  6.28) 

So  (7.4) 

in  which  the  right  side  is  a function  of  x and  y alone.^^  Further  in  equation  7.1, 
z is  a function  of  infiltration  time  alone  and  hence,  for  any  given  x and  t,  is 
dependent  on  the  past  history  of  advance  only.  Indeed,  then 


^ = i(^)  (7.5) 

in  which  i is  the  infiltration  rate,  a function  of  t,  the  infiltration  time.  Thus, 
equation  7.1  can  be  viewed  as  a partial  differential  equation  in  a single 
unknown,  say  y,  in  terms  of  two  independent  variables  x and  t. 


c 


dX 


+ 


at 


+ i = 0 


(7.6) 


in  which 


P _ dq(y) 
' dy 


(7.7) 


or  with  equation  7.4  in  force. 


c = ( y3'2  + y Co  y''2  ) 


(7.8) 


Finally,  with  the  Chezy  Co  expressed  in  terms  of  the  Manning  n (eq.  6.31) 


c 


i A 

3 n 


(7.9) 


As  for  the  zero-inertia  model,  the  program  for  kinematic-wave  solutions  works 
with  dimensionless  variables.  The  techniques  for  reducing  variables  and 
equations  to  dimensionless  form  are  the  same  as  described  in  section  6.8. 
The  equations  in  this  chapter  have  the  same  form  regardless  of  whether  the 
symbols  stand  for  dimensioned  variables  or  their  dimensionless  counterparts 
and  can  be  viewed  either  way. 


25.  . .y  alone,  in  the  event  of  a constant  bottom  slope  ancf 
roughness. 
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7.2  Reduction  to  a Pair 
of  Ordinary  Differential 
Equations 


Equations  like  equation  7.6  describe  a kind  of  wave  motion  called  kinematic, 
in  the  nomenclature  of  Lighthill  and  Whitham  (3),  pioneers  in  the  area.  The 
name  stems  from  the  fact  that  all  the  variables  in  the  governing  equation  (7.6) 
are  kinematic  in  nature;  that  is,  their  dimensions  involve  length  and  time  only. 
In  a sense,  the  nomenclature  is  unfortunate,  because  equation  7.6,  besides 
expressing  the  equation  of  continuity,  which  is  certainly  a kinematic  relation, 
has  also  implied  equation  7.3,  a dynamic  relation,  stating  equality  of  weight 
and  drag  force  components. 

This  equality,  describing  uniform  flow  at  normal  depth  at  every  point  in  the 
stream  is  met,  approximately,  whenever  the  depth  gradient  is  small  compared 
with  the  bottom  slope.  This  will  be  true  in  the  bulk  of  the  surface  stream 
when  the  bottom  slope  is  relatively  large.  The  method  fails  on  the  spot,  as 
soon  as  a zero  or  adverse  bottom  slope  is  encountered. 

At  the  very  front  of  an  advancing  stream,  and  also  at  the  brink  of  a free  over- 
fall,  the  depth  gradient  becomes  very  large,  and  there  equation  7.3  is  certainly 
not  satisfied.  Furthermore,  behind  an  end  block,  where  the  ponded  water  is 
essentially  stagnant,  even  on  a steep  slope,  the  depth  gradient  and  bottom 
slope  are  essentially  equal  in  magnitude.  On  the  other  hand,  the  condition  of 
equation  7.3  is  very  well  met  in  the  trailing  region  of  the  surface  stream  in 
recession,  a region  that  is  generally  difficult  to  compute  on  steep  slopes  with 
the  zero-inertia  model,  because  of  the  low  depths  there. 

In  general,  limits  to  the  reasonable  application  of  kinematic-wave  theory  to 
border  irrigation  can  be  found  by  comparing  the  results  thereof  with  those 
stemming  from  the  zero-inertia  model.  Within  those  limits  the  zones  of  large 
negative  depth  gradient  are  effectively  modeled  by  depth  discontinuities. 

Also,  behind  a blocked  border  end,  the  water  can  be  assumed  a stagnant 
pool,  having  no  effect  on  the  stream  plunging  in.  The  stream  calculations, 
performed  with  the  kinematic-wave  assumption,  equation  7.3,  in  force  are 
then  simply  terminated,  in  effect,  at  the  point  in  the  border  marking  the 
upstream  boundary  of  the  pond;  the  volume  of  the  pond  is  determined  by  sub- 
tracting the  infiltrated  volume  seeping  out  of  the  pond  from  the  stream 
volume  entering. 

These  techniques  are  detailed  in  the  rest  of  this  chapter,  following  an  exposi- 
tion in  the  next  sections  of  the  basic  approach  in  solving  equation  7.6. 

By  an  appropriate  change  of  independent  variables,  equation  7.6  can  be 
transformed  to  a pair  of  ordinary  differential  equations,  which  can  then  be 
solved  numerically.  The  first  step  consists  in  introducing  an  orthogonal  s-n 
system  oriented,  at  each  point  in  the  x-t  plane,  at  some  angle  a (as  yet 
undefined)  to  the  x-t  coordinates.  In  general,  a will  vary  from  point  to  point  in 
the  plane,  as  in  figure  7.1,  so  that  the  s and  n coordinates  form  two  mutually 
orthogonal  families  of  curves. 


76 


X 

Figure  7.1.  — The  s-n  coordinate  system. 


By  the  chain  rule  of  differentiation,  equation  7.6  is  transformed  to 


ay  as  ay  an  ay  as  ay  an  . ^ 

c — + c — + — + — +1  = 0 


as  ax 


an  ax 


as  at 


an  at 


(7.10) 


With 


as  an 

= cos  a 


ax 


ax 


3s  . an  ^ „ /-7  H 

- sin  a = sin  a = cos  a (7.11) 


at 


at 


equation  10  can  be  written 


(c  cos  a + sin  a)  — (c  sin  a - cos  a)  +1=0  (7.12) 

as  an 


If  now  a is  chosen  at  each  point  in  the  x-t  plane,  in  accord  with  the  value  of  c 
that  in  fact  exists  there,  such  that 


cot  a = c (7.13) 

the  coefficient  of  ay/an  in  equation  7.12  vanishes,  and  this  equation  (since 
c = cos  a/sin  a,  and  cos^  a.  + sin^  a = 1)  can  be  written  as  the  ordinary 
differential  equation. 


— ' ^ + i = 0 (7.14) 

sin  a ds 
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But 


1 dy  _ 
sin  a ds  dt 


(7.15) 


as  can  be  seen  from  figure  7.2,  and  so  equation  7.14  takes  the  simple  form 


dt 


(7.16) 


valid  along  a so-called  characteristic  curve,  whose  equation  is  given  by 


(7.17) 


Equation  7.17  follows  from  equation  7.13,  when  the  relations  depicted  in 
figure  7.2  are  taken  into  account.  Equations  7.16  and  7.17  represent  the 
ordinary  differential  equations  equivalent  to  the  single  partial  differential 
equation  7.6. 


The  x-t  plane  can  be  thought  of  as  blanketed  with  a family  of  these  character- 
istic curves  (the  s-coordinate  curves),  for  which  equation  7.17  expresses  the 
slope  at  any  point  in  terms  of  the  depth  there.  Along  any  such  curve,  the 
depth  decreases  with  time  as  the  result  of  infiltration,  in  accord  with  equa- 
tion 7.16.  The  general  appearance  of  these  curves  will  be  depicted  in  section 
7.4,  once  an  important  complication,  intersection  of  characteristic  curves,  has 
been  dealt  with. 
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7.3  Intersection  of 
Characteristics:  The 
Kinematic  Shock 


The  formation  of  kinematic-shock  discontinuities  will  be  studied  first  for  a 
special  case— the  gradual  increase  of  a preexisting  uniform  flow  in  a channel 
with  zero  infiltration  and  constant  roughness  and  bottom  slope.  The  results 
will  then  be  viewed  in  the  light  of  conditions  that,  in  fact,  exist  in  a border- 
irrigation  flow. 

Figure  7.3  depicts  the  characteristic  curves,  or  simply  characteristics,  per- 
taining to  an  initially  uniform  flow  in  a reach  of  impervious  channel,  steep 
enough  for  the  wave  motion  therein  to  be  described  by  kinematic-wave 
theory.  At  time  zero,  let  the  discharge  gradually  increase.  This  simplified  case 
retains  the  essential  features  necessary  to  extend  the  results  to  the  case  of 
flow  with  infiltration. 

In  the  given  special  case,  the  characteristics  are  all  straight  lines,  each  bear- 
ing a single,  constant  value  of  depth  and  discharge;  in  the  region  of  uniform 
flow  the  characteristics  are  parallel.  Their  slope,  on  t = 0,  is  given  by 
equations  7.17  and  7.18.  At  x = 0,  the  increasing  discharge  is  accompanied 
by  an  inverse  slope  increasing  with  time.  Consequently,  the  characteristics 
intersect.  The  significance  of  the  intersection  will  be  dealt  with  in  due 
course.  Of  immediate  concern  are  the  conditions  governing  the  time  and 
place  of  intersection  of  a neighboring  pair  of  characteristics,  as  depicted  in 
figure  7.4. 


Figure  7.3.  — Formation  and  propagation  of  kinematic 
shocks. 
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Figure  7.4.  — Intersection  of  a neighboring  pair  of 
characteristics. 


Of  course,  the  smaller  is  5t,  the  more  nearly  are  the  characteristics, 
emanating  from  x = 0 at  t = 0 and  t = 6t,  parallel.  At  the  same  time,  the 
smaller  is  their  initial  separation,  so  that  intersection  is  more  easily  achieved. 
It  turns  out  that  as  6t  approaches  zero,  limiting  values  of  the  coordinates  (Xg, 
tg)  of  intersection  are  approached.  Indeed,  with  c and  dc/dt  the  inverse  slope 
and  its  rate  of  increase,  respectively,  at  x = 0 at  time  zero,  then  c + (dc/dt)6t 
is  the  inverse  slope  of  the  characteristic  emanating  from  x = 0 at  t = 6t.  The 
triangles  constructed  in  figure  7.4  show  that 

Xg  dC  Xg 

c = — and  c + — - 6t  = — (7.18) 

tg  dt  tg  - 6t 

With  Xg  eliminated  from  equations  7.18,  and  6t  made  smaller  and  smaller,  it  is 
evident  that  in  the  limit. 


dc/dt 


(7.19) 


dc/dt 


(7.20) 
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From  these  results,  it  is  evident  that  an  abrupt  increase  in  discharge  (and 
hence  c)  at  the  left  boundary  will  result  in  immediate  formation  of  a shock. 
Further,  it  is  clear  that  as  the  initial  depth  (and  c)  decreases,  the  intersection 
occurs  ever  earlier.  In  the  limit,  even  a gradual  increase  in  discharge,  from 
zero,  will  lead  to  immediate  shock  formation. 

Since  the  only  effect  of  infiltration  and  variable  slope  and  roughness  is  to 
give  the  characteristics  some  curvature,  these  qualitative  results  are 
unchanged^®  for  an  irrigation  stream  in  a pervious  border.  Thus,  the  motion  of 
the  surge  onto  the  field  commences,  from  the  start,  as  a kinematic  shock. 

The  significance  of  this  is  discussed  next. 

Again,  the  argument  is  presented  for  the  simplified  case— uniform,  imper- 
vious channel.  The  qualitative  conclusions  are  not  altered  for  the  flow  in 
border  irrigation.  It  has  already  been  pointed  out  that  each  characteristic 
carries  its  own  constant  value  of  depth.  Thus,  as  two  characteristics 
approach  each  other,  the  longitudinal  distance  between  the  two  different 
depth  values  decreases.  As  the  intersection  occurs,  the  heretofore  steepening 
profile  becomes  actually  vertical,  that  is,  a discontinuity.  This  abrupt  step 
change  in  depth  constitutes  the  profile  of  the  kinematic  shock.^^ 

A shock-discontinuity  propagates  at  a velocity  different  from  the  velocity  of 
either  of  the  kinematic  waves  associated  with  the  depths  on  its  high  and  low 
sides,  and  turns  out  to  lie  in  between  these  values.  Elementary  mass- 
conservation  applied  to  a thin  slice  of  space  through  which  the  shock  is  mov- 
ing, as  in  figure  7.5,  shows  that  shock-propagation  velocity  Wg  is  the  same  as 
velocity  Va  of  the  water  just  behind  the  wave  front. 

Indeed,  let  be  a plane  fixed  in  space  through  which  water  is  flowing  and 
which  is  located  a very  small  distance  6x  behind  the  shock  front.  Let  Xg  be  a 
moving  plane,  fixed  in  the  shock  front  of  height  yg.  At  the  given  instant,  the 
volume  of  water  between  x^  and  Xg  is  growing  at  the  rate  Wg  ya  plus  small 
quantities  associated  with  the  infiltrated  volume  and  the  gradually  changing 
depth  behind  the  shock.  The  volumetric  rate  of  flow  into  the  region  between 
Xi  and  Xa  is  Vya,  plus  a small  quantity  associated  with  the  depth  gradient  and 


26As  very  small  Xg  are  considered,  any  curvature  short  of  in- 
finite will  not  be  apparent  and  infinite  curvature  cannot  arise 
even  at  the  wave  front  where  the  infiltration  is  infinite, 
because  the  surface  depth  never  drops  to  zero,  for  otherwise 
the  velocity  must  also. 

27.  . .that  is,  as  viewed  in  this  study,  in  which  Sq  = Sf  is 
strictly  observed.  It  is  possible  (see,  for  example,  Henderson 
(2))  to  obtain  a more  realistic  shock-front  profile  by  relaxing 
this  condition  in  the  vicinity  of  the  front,  allowing  instead 
both  a depth  gradient  and  local  and  convective  accelera- 
tions, assuming  only  a profile  slowly  changing  in  shape. 
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Figure  7.5.  — Shock  propagation  speed. 


7.4  Pattern  of 
Kinematic  Waves  in 
Border  Irrigation 


6x  between  Xi  and  Xg.  But  the  volumetric  inflow  rate  at  Xi  must  equal  the  rate 
of  growth  of  volume  between  x^  and  Xg.  Further,  as  6x  is  considered  smaller 
and  smaller,  the  small  correction  terms  become  negligible.  Thus, 


dXg 

dt 


(7.21) 


in  which  Vg,  qg,  and  Vg  are  the  values  of  V,  q,  and  y,  respectively,  when  6x  — 0. 


The  front  of  the  irrigation  stream,  then,  moves  out  from  the  origin  of  the  x-t 
plane  as  a depth  discontinuity.  The  shock  trajectory  is  intersected  by  suc- 
cessive kinematic  waves  emanating  from  the  left  boundary.  At  each  inter- 
section, the  depth  on  the  high  side  of  the  shock  is  that  carried  by  the 
kinematic  wave;  the  shock-propagation  speed  leaving  the  intersection  reflects 
the  new  depth  and  discharge  on  its  high  side.  A complete  pattern  of  trajec- 
tories is  shown  in  figure  7.6.  Visible  are  the  intersections  of  the  advancing 
shock  trajectory  with  a succession  of  kinematic-wave  trajectories,  each 
emanating  from  the  upstream  field  boundary  at  an  ever  greater  value  of  time. 
The  arrival  of  waves  at  the  right  boundary  subsequent  to  the  arrival  of  the 
stream  front  there  yields  runoff  from  field  end.  Upon  cutoff,  a fan  of 
kinematic-wave  trajectories  radiates  outward  from  the  cutoff  point  (x  = 0, 
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t = tco).  Some  reach  field  end  to  contribute  to  runoff.  The  rest  terminate  (with 
zero  depth)  in  the  interior  of  the  field;  their  end  points  constitute  the  reces- 
sion curve. 

The  border  utilized  in  the  construction  of  figure  7.6  was  uniform  in  slope  and 
roughness,  and  the  inflow  rate  was  held  constant  until  cutoff.  The  kinematic- 
wave  trajectories  curve  upwards  slightly  because  of  the  effect  of  infiltration 
along  their  lengths  (eqs.  7.16,  7.17,  7.8);  the  advance  trajectory  also  curves 
upward,  reflecting  the  decreasing  depth  on  the  upstream  side  of  the  shock. 

When  the  shock  reaches  field  end,  advance  is  over.  Successive  kinematic-wave 
trajectories  each  bring  a value  of  discharge  to  field  end;  this  constitutes  the 
runoff  (or  rate  of  ponding  if  a dike  is  assumed  at  field  end).  In  the  event  that 
high-side  shock  depth  reduces  to  zero  before  reaching  field  end,  advance 
halts  spontaneously  and,  usually,  front-end  recession  follows. 

Recession  from  the  upstream  end  commences  immediately  after  cutoff,  as 
discharge  and  depth  fall  to  zero  simultaneously  in  the  normal-depth  kinematic- 
wave  approximation.  However,  during  the  instant  that  the  depth  is  falling, 
kinematic  waves  are  still  being  generated,  each  at  a successively  lower  value 
of  inverse  slope.  Each  of  these  waves  is  subject  to  the  conditions  of  equation 
7.16.  Some  of  them,  those  with  the  larger  starting  depths,  may  yet  get  to  field 
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Figure  7.6.  — Complete  network  of  characteristic  curves. 


83 


end  and  provide  some  runoff.  Those  with  smaller  starting  depths  end  spon- 
taneously, somewhere  in  the  border,  as  the  depth  upon  them  reduces  to  zero 
because  of  infiltration.  These  end  points  mark  the  trajectory  of  the  recession 
curve. 


7.5  Program 
Organization 


7.5.1  Definition  of  Mesh  Density 


The  denseness  of  the  network  of  characteristic  curves  originating  prior  to 
cutoff  is  governed  by  the  time  increments  6t  separating  the  characteristics 
emanating  from  the  upstream  boundary.  This  is  chosen  on  the  basis  of  the 
input  variable  Nstd,  entered  in  line  10  (sec.  3.10), 


(7.22) 


in  which  Wg  and  c are  computed  for  the  initial  inflow  rate.  This  makes  the 
number  of  points  defining  the  advance  curve  somewhat  greater  than  Nstd-^® 

7.5.2  Construction  of  First  Increment  of  Advance 

The  first  step  in  constructing  the  network  of  characteristics  is  accomplished 
by  a call  to  the  subroutine  KWCDA,  which  is  described  in  detail  in  section  7.6. 
In  brief,  it  computes  the  intersection  of  a segment  CA  of  a kinematic-wave 
trajectory  with  a segment  DA  of  the  shock  trajectory  (advance  curve).  In  this 
case,  C and  D are  both  on  x = 0,  (see  fig.  7.7)  and  are  separated  by  the  time 
increment  5t. 

The  subroutine  returns  the  coordinates  (Xg,  tj  of  point  A,  and  the  depth  yg, 
discharge  Pg,  celerity  c,  and  advance  speed  Wg  there. 

7.5.3  Construction  of  the  Family  of  Characteristics  Before  Cutoff 

Subsequently,  the  program  constructs  each  successive  kinematic-wave  trajec- 
tory in  a series  of  distance  steps.  The  node  points  have  the  same  x-values  as 
the  node  points  on  the  advance  curve  below.  When  a given  kinematic-wave 
trajectory  has  been  constructed  out  to  the  current  limits  of  the  advance 
curve,  the  subroutine  KWCDA  is  called  to  construct  both  the  final  segment  of 
the  kinematic-wave  trajectory  and  the  final  segment  of  the  advance  curve. 
Their  intersection  provides  a new  node  location  for  the  mesh. 

It  is  presumed  that  as  the  kinematic  wave  trajectories  are  constructed,  one  by 
one,  separated  at  their  upstream  end  by  the  time  increment  6t  (eq.  7.22),  even- 
tually one  will  be  encountered  for  which  the  last  segment  straddles  field  end. 


2®With  zero  infiltration  and  constant  bottom  slope  and 
roughness,  it  would  be  just  Nsjd  -i-  1. 
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Figure  7.7.  — First  step  in  construction  of  network  of 
characteristics. 


The  time  of  arrival  of  the  stream  at  field  end  is  determined  by  interpolation 
along  the  advance  curve.  The  value  of  runoff  (or  rate  of  increase  of  ponded 
volume),  that  is,  the  discharges  at  x = L on  both  the  kinematic-wave  trajec- 
tory and  the  advance  trajectory  are  found  by  equation  7.4  after  the  depths 
have  been  found  there  by  interpolation. 

Henceforth,  at  least  until  cutoff,  field  end  constitutes  the  x-location  of  the 
downstream-most  node  in  the  network  of  characteristic  curves.  The  net  is 
built  up  through  repeated  use  of  KWBC  on  each  characteristic;  the 
downstream-most  value  of  discharge  on  each  constitutes  the  rate  of  runoff  or 
increase  in  ponded  volume.  Each  characteristic  is  identified  by  its  index  i 
(which  should  not  exceed  100  to  stay  within  array  dimensions). 
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7.5.4  Post  Cutoff  Family  of  Characteristics 


After  inflow  cutoff,  a number  N^ec  of  kinematic-wave  trajectories  is  chosen 
sufficiently  large  to  define  the  fan  emanating  from  x = 0,  t = tgo.  In  general, 
those  members  of  the  fan  bearing  the  largest  discharge  values,  initially,  will 
reach  the  right  boundary.  The  rest  will  turn  upward  and  end  spontaneously, 
within  the  border  length,  as  the  discharge  thereon  reduces  to  zero.  The 
number  of  characteristics  in  the  fan  is  given  by  the  formula 


N,ec  = 2 + —^ 


(7.23) 


6t  . c 


in  which  c^  is  the  celerity  associated  with  the  discharge  before  cutoff.  This 
number  has  proved  to  provide  a network  after  cutoff  comparable  in  density  to 
the  portion  before  cutoff. 

To  ensure  a reasonably  uniform  distribution  of  members  within  the  fan,  the 
following  formula  provides  the  starting  value  of  discharge  q(1)  for  each 
member. 


(7.24) 


in  which  is,  for  the  moment,  the  precutoff  value  of  inflow,  and  j^ec  is  an 
index  for  characteristics  within  the  fan. 

The  trajectory  of  each  member  of  the  fan  is  computed  by  repeated  use  of 
KWBC.  Each  characteristic,  as  it  reaches  the  last  node,  at  field  end,  provides 
a value  of  discharge,  for  runoff  or  ponded-volume  increase.  Sooner  or  later, 
however,  a characteristic  will  be  encountered  that  fails  to  reach  field  end  (a 
negative  depth  will  be  computed). 

At  this  time  a subroutine,  KWEND,  is  introduced.  KWEND  is  designed  to 
compute  a kinematic-wave  trajectory  all  the  way  to  its  end  (y  = 0),  given  the 
location  and  discharge  at  its  starting  point.  This  subroutine,  described  in 
detail  in  section  7.8,  chooses  its  own  node  locations,  independent  of  the 
node  distribution  on  the  advance  curve.  This  proves  necessary  because  of  the 
great  curvature  and  small  inverse  slopes  encountered  on  the  characteristics 
as  their  trailing  ends  (y  = 0)  are  approached  (see  fig.  7.6). 
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First,  KWEND  is  used  repeatedly,  in  an  iterative  scheme  to  find  that  value  of 
discharge  q(1)  at  the  starting  point  (x  = 0,  t = t^o)  of  a characteristic  such 
that  its  end  point,  y = 0,  falls  right  on  field  end,  x = L.  A discharge  greater 
than  this  critical  value  yields  a characteristic  with  its  end  value,  Xgnd,  greater 
than  L,  while  one  less  than  this  value,  results  in  a characteristic  ending  at  an 
X less  than  In  the  vicinity  of  the  correct  value  of  q(1),  the  relation  between 
q(1)  and  Xgnd  is  essentially  linear;  this  provides  relatively  quick  convergence  of 
the  iterative  scheme.  The  time  coordinate  of  the  end  point  of  the  critical 
characteristic  marks  the  end  of  the  runoff  function  qro(t),  and  also  the  last 
point  on  the  recession  function  x,ec(t). 

At  this  point  in  the  program,  is  redefined 

Nrec  = Nstd  (7.25) 

and  jrec  is  reset  to  zero.  At  the  same  time  qi  in  equation  7.24  is  set  to  the 
critical  value  q(1)  just  found. 

Subsequently,  equation  7.24  is  used  to  yield  starting  discharge  values  for 
each  of  N,ec  characteristics  (as  j,ec  is  incremented  by  one  for  each  successive 
characteristic),  the  end  points  of  which  will  define  the  recession  curve  at  x 
values  interior  to  0 and  L.  The  subroutine  KWEND  constructs  each  such 
characteristic  and,  in  particular,  the  x and  t coordinates  of  its  end  point. 

When  all  the  characteristics  have  been  computed  and  all  end  points  have 
been  joined  to  form  the  recession  function,  the  program  branches,  depending 
on  the  given  downstream  boundary  condition. 

With  DBG  = 1 entered  as  input  data  (sec.  3.5),  to  signify  free  flow  off  border 
end,  the  computation  of  runoff  and  recession  are  accepted,  and  the  program 
proceeds  to  post-irrigation  analysis  (sec.  5.5).  With  DBG  = 2 (ponding),  the 
program  computes  the  elevation  of  a pool  at  field  end  whose  volume  equals 
that  computed  as  runoff.  The  upstream-most  extent  of  that  pool  marks  the 
downstream  limit  of  validity  of  the  recession  curve  as  computed  earlier.  The 
rest  of  the  recession  curve  is  calculated  as  resulting  from  infiltration  of  a still 
pond.  When  this  computation  is  completed,  the  program  proceeds  to  post- 
irrigation analysis  (sec.  5.5). 


29a  value  of  IDIAG  > 4 will  result  in  a printout  describing 
each  trial  characteristic;  IPLOTC  > 2 will  result  in  a plot  of 
each  trial  trajectory. 
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7.5.5  The  Hybrid  Model:  Recession,  Only,  Computed  by  Kinematic-Wave  Theory 

With  the  option  SOLMOD  = 5 in  force,  the  flow  is  calculated  by  zero-inertia 
theory  until  the  stream  has  reached  field  end  and  also  until  recession  has 
started  from  the  upstream  end.  Then,  the  surface-depth  profile  computed  for 
that  time,  called  t^wsu,  is  modified  for  transfer  to  kinematic-wave  theory,  and 
the  rest  of  the  irrigation  is  computed  on  that  basis. 

The  depth  values  in  the  profile  are  adjusted  by  averaging  neighboring  values 
to  assure  a monotonically  increasing  depth  with  distance.  In  the  front,  where 
the  depth  gradient  is  basically  negative,  the  depth  is  held  constant.  Surface- 
water  volume  is  preserved  by  shortening  the  stream  as  necessary. 

A kinematic-wave  trajectory  is  then  assumed  to  emanate  from  each  node  on 
the  time  line  t = tk^^su,  with  the  starting  conditions  derived  from  the  depth  per- 
taining to  that  node.  The  subroutine  KWEND  is  called  to  construct  each 
characteristic.  If  any  characteristic  proves  to  cross  field  end  instead  of  ter- 
minating within  the  border,  the  time  of  crossing  and  the  discharge  at  that 
time  are  found  by  interpolation  to  contribute  to  the  runoff  function. 

As  described  in  section  7.5.4,  the  recession  curve  is  found  by  joining  the  ends 
of  the  kinematic-wave  trajectories  which  terminate  within  the  border.  The 
passage  of  the  recession  curve  through  field  end  is  found  by  interpolation 
between  the  two  termination  points  closest  to  field  end,  one  on  each  side 
(fig.  7.8). 

7.6  Subroutine  KWCDA  This  subroutine  is  called  every  time  that  the  intersection  of  a segment  of  a 

characteristic  and  a segment  of  the  advance  trajectory  must  be  computed. 
Thus,  it  is  called'to  compute  the  last  segment  of  every  characteristic  during 
the  advance  phase.  It  performs  a numerical  integration  of  equations  7.16, 

7.17,  and  7.21  over  segments  of  the  characteristic  and  advance  trajectories 
pictured  in  figure  7.9. 

The  integration  of  equation  7.16, 


is  not  trivial,  because 
the  problem.  Indeed, 


Ya  - Yc  = - 


(7.26) 


00,  as  t — t^.  A change  of  variables,  however,  solves 


fA  ~ fc 

TC 


Zc 


(7.27) 
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Figure  7.8.  — Characteristic  curves  in  hybrid  model. 


X 


Figure  7.9.  — Intersection  of  characteristic  and  advance 
curves. 
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Thus,  a numerical  approximation  to  equation  7.26  is  given  by 


Yc  - Va  = ^ (tA  - tc)  (7.28) 

TC 

Equation  7.28,  which  has  and  tA  as  unknowns,  is  solved  simultaneously 
with  the  following  approximations  to  equations  7.17  and  7.21,  respectively, 

Xa  — Xc  _ Cc  + Ca 

^A  ~ fc  2 

Xa  - Xp  _ Wp  + Wa 

^A  ~ fo  2 


(7.29) 

(7.30) 


which  introduce  the  additional  unknown  Xa.  There  are  no  further  unknowns, 
because  Ca  and  Wa  are  uniquely  tied  to  yA  through  equations  7.8,  7.4  and  the 
last  of  equations  7.21. 

If  the  linear  mode  of  solution  is  in  effect  (that  is,  LINMOD  = 1— sec.  3.9),  the 
segments  in  figure  7.9  are  assumed  straight  lines,  so  that  Cq  = Ca,  Wp  = Wa, 
and  a direct  solution  is  possible.  With  the  nonlinear  option  (LINMOD  =2). 
the  equations  are  solved  iteratively.  If  the  iterations  fail  to  converge  in  20 
steps,  the  entire  set  of  iterations  is  repeated  with  the  results  of  each  step 
printed  to  aid  in  a diagnosis  of  the  problem. 


7.7  Subroutine  KWBC  The  aim  of  this  subroutine  is  to  construct  a segment  of  kinematic-wave  tra- 
jectory. Given  the  conditions  at  its  start,  at  point  B in  figure  7.10,  namely,  Xp, 
^B.  Yb.  tq,  Zp,  qp,  Cp,  the  x-coordinate  of  its  end  point  Xq,  and  the  time  coor- 
dinate tp  of  the  advance  curve  at  the  point  Xp  = Xp,  KWBC  determines  the 
remaining  conditions  at  the  end  point,  tp  , Yc,  rp,  Zp,  qp,  Cp. 

This  is  achieved  by  a numerical  integration  of  the  equations  7.16  and  7.17, 
that  is,  by  simultaneous  solution  of  algebraic  approximations  to  their 
integrals,  namely. 


Yc  ~ Yb  _ Zp  — Zp 
fc  ~ fs  Tc  ~ T% 


(7.31) 


Xp  ~ Xp  Cp  -I-  Cp 
~ ^B  2 


(7.32) 


Equation  7.31  is  derived  from  equation  7.16  in  a manner  essentially  the  same 
as  that  leading  to  equation  7.28.  Equation  7.32  is  the  same  as  equation  7.29 
in  meaning.  With  rp  = tp  - tp,  and  Zp  and  Cp  unique  functions  of  Tp  and  yp, 
respectively,  equations  7.31  and  7.32  are  solved  simultaneously  for  the 
unknowns  tp,  yp. 
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7.8  Subroutine  KWEND 


X 


Figure  7.10.  — Construction  of  a segment  of  the  kinematic- 
wave  trajectory. 


In  the  linear  mode  (LINMOD  = 1— sec.  3.9),  it  is  assumed  that  BC  is  a 
straight  line,  and  hence  that  Cg  = Cq.  This  permits  direct  solution  of  equation 
7.32  for  tc.  followed  by  another  direct  solution,  of  equation  7.31,  for  Vc.  In  the 
nonlinear  mode  (LINMOD  = 2),  BC  is  considered  a segment  of  the  second- 
degree  curve  described  by  equation  7.32  with  Cg  ^ Cq,  and  an  iterative 
Newton-Raphson  solution  is  required.  The  linear  solution  then  constitutes  the 
first  guess.  If  the  iterations  fail  to  converge  within  20  cycles,  all  iterations  are 
repeated  with  the  results  of  each  cycle  printed  as  a diagnostic. 

This  subroutine  requires  given  x and  t coordinates  and  a given  discharge  at 
that  place  and  time.^°  The  subroutine  constructs  a kinematic-wave  trajectory 
(characteristic)  emanating  from  the  given  point  and  bearing  the  given  value  of 
discharge.  The  trajectory  is  constructed  to  its  end,  where  the  discharge  and 
depth  have  fallen  to  zero  as  the  result  of  infiltration.  Knowledge  of  the 
advance  curve  and  infiltration  function  is  a requirement  for  the  use  of 
KWEND. 


3°.  . .and,  hence,  depth,  celerity,  infiltration  depth,  and  so 
forth. 
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The  characteristic  is  constructed  in  time  steps,  by  numerical  integration  of 
equations  7.16  and  7.17  over  each  step.  The  size  of  the  step  is  given  before 
the  successful  construction  of  each  segment.  At  first,  the  step  is  based  on 
the  assumption  that  the  total  length  of  characteristic  will  be  composed  of 
N = Nstd^^  segments.  The  segments  are  numbered  in  order,  with  the  index 
k = 2 referring  to  the  first  segment.^^  The  label  B is  given  to  the  start  of  the 
segment  and  C to  the  end  of  the  segment. 

At  the  start  of  the  calculation  for  each  segment,  an  estimate  is  made  of  the 
time  tend  at  which  the  entire  trajectory  will  end,  namely, 

tend  = to  + rfZend)  (7.33) 

In  equation  7.33,  tp  is  the  time  at  which  the  advance  curve  arrived  at  the  point 
in  the  field  Xg,  at  which  the  current  segment  starts,  and  rfZend)  is  the  time  it 
would  take,  with  the  given  infiltration  function,  to  infiltrate  a depth 
Zend  = Zb  + ye-  Clearly,  the  more  nearly  vertical  is  the  remainder  of  the 
characteristic  leaving  point  B,  and  the  more  nearly  horizontal  the  advance 
curve,  the  more  nearly  correct  will  this  estimate  be. 

The  projected  size  of  the  time  step  encompassed  by  the  segment  under  con- 
struction is  given  by  the  formula 


6t 


fend  fs 

N - k -I-  2 


(7.34) 


and  then. 


tc  — fe  ^f  (7 .35) 

An  attempt  is  made  to  construct  the  segment  by  simultaneous  solution  of 
numerical  approximations  to  the  integrals  of  equations  7.16  and  7.17,  namely. 


Vc  ~ Vb  _ Zc  — Zb 
fc  ~ fe  Tc  ~ t’b 

which  is  the  same  as  equation  7.31,  and 

Xq  ~ ^B  _ Cb  -k  Cc 

fc  ~ fB  2 


(7.36) 


(7.37) 


3'' . . .entered  as  data  input  (see  sec.  3.10). 

32k  is  also  the  index  of  the  node  points  on  the  characteristic. 
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a restatement  of  equation  7.32.  Equation  7.36  suffers  from  a potential  prob- 
lem. Quite  possibly,  at  some  point  in  its  length,  the  kinematic-wave  trajectory 
will  parallel  the  advance  curve  below.  Then  the  denominator  of  the  expres- 
sion in  the  right  side  will  vanish.  The  program  is  prepared  for  this  even- 
tuality— the  denominator  is  monitored  and  when  its  absolute  value  falls  below 
0.000001,  equation  7.36  is  replaced  by 


Vc  - Yb 
t - t “ “ 

with  i computed  at  the  current  trial  location  of  C. 


(7.38) 


In  the  linear  mode  of  solution  (LINMOD  = 1— see  sec.  3.9),  equation  7.37  is 
solved  first,  for  Xc,  in  the  assumption  that  Cc  = Cq.  With  tp  the  time  the 
advance  curve  crosses  Xq,  Zc  = zfrj  is  found  for 


rc  = tc  - to  (7.39) 

This  allows  solution  of  equation  7.36  (or  7.38)  for  yc. 

With  LINMOD  = 2 (the  nonlinear  solution  mode),  the  linearized  solution 
forms  the  first  guess  in  an  iterative  solution  of  equations  7.36  (or  7.38)  and 
7.37  for  Xc  and  yc.  If  20  iterations  prove  insufficient  for  convergence,  all  20  are 
repeated  and  the  results  printed  for  diagnosis  of  the  problem.  In  the  event 
that  yc  turns  negative,  6t  is  cut  in  half  and  the  sequence  of  instructions  com- 
mencing with  equation  7.35  is  repeated. 

When  all  N segments  have  been  completed,  a final  end  point,  upon  which 
y = 0,  is  computed  in  accordance  with  the  equations 

x.„d  = Xc  + (7.40) 

'end 

Ua  = tc  + — (7.41) 

'end 


in  which  C represents  the  end  of  the  last  segment  computed.  Equation  40  is 
an  approximation  to  the  integral  of  equation  7.17,  for  with  equations  7.7  and 
7.16  in  force, 


^end 


'end  dq  dy/dt 
'c  dy  i 


dt 


Equation  7.41  is,  of  course,  an  approximation  to  the  integrall  of  equation  7.16. 
In  both  equations  7.40  and  7.41,  iend 's  an  approximation  to  i. 
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Appendix:  Example  Input/Output 


Printout  from  a sample  run  is  presented. 

Input  keyboard  entries  are  marked  by  an  arrow.  Immediately  after  the  entry  is 
completed,  the  program  echoes  back  the  information  received.  This  fulfills 
two  purposes:  (1)  It  verifies  that  the  intended  data  were  received  and  (2)  if  the 
program  is  used  in  batch  mode  with  punched-card  input,  the  echoed  data 
following  the  printed  prompts  provide  a record  of  the  data  input. 
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ttttUtttt>tttltt>lttlttttttttBORDER-IRRIUTION-FLOU  PROGRAHtttlIttttttmtttttmtmtm 


08  JUN  83  18.17.30 

aAPSED  CP  TIME  (SECONDS)  = 35.0A5 
CU-S  USED  = 115.276 
CU-S  REWlINING  = 584.724 


INTERACTIVE  USERS  --  POSITION  PAPER  FOR  PROMPTS  TICN  ENTER  (LINE  1)  1 FOR  A FRESH  STWIT  -- 

2 TO  CtMIGE  DESIGN  PMtAiCTERS  - ZREQiOINtTCO  — 

3 TO  0MN6E  SOLUTION  MOOES  - SOLMODfLINMODfOTHOD*  ISUPZA>ZADMOD  — 

4 TO  CHANGE  SOLUTION  PARAMETERS  - NSTDiRDXiDTSTDfTMAXf JMAX  — 

5 TO  CHANGE  LlVil  OF  DIAGNOSTICS  AND/OR  PLOTTING  FLAGS  — 

0  TO  STOP  — 


1  

1 

(LINE  2) IDENTIFY  THIS  RUN  

EXAMPLE  11  

EXAMPLE  41 

ENTER  (LINE  3)  --  INPMOD  DMLMOD  — 

2,1  

2  1 

ENTER  (LINE  4)  SOIL  AND  CROP  HYDRAULIC  PROPERTIES  --  RUFMOD  RUF  AN  INFMOD  K A B C — 

2t  0.04f  0.  1.  1>  0.5.  0.  0.  -« 

2 .04  0.  1 1.  .5  0.  0. 

ENTER  (LINE  5)  FIELD  GEOMETRY  --  L DBC  SOMOD  — 

1200.  1.  1 

1200.  1 1 

ENTER  (LINE  6)  AVERAGE  BOTTOM  SLOPE  - SOAVG  — 

0.0005  « 

.0005 

ENTER  (LINE  7)  MANAGEMENT  PARAMETERS  --  ZREQ  0 TCO  ~ 

1.5.  0.2.  40  

1.5  .2  40. 

ENTER  (LINE  9)  SOLUTION  PARA«TERS  --  SOLMOD  LINMOD  DTMOD  ISUPZA  ZADMOD  — 

2.  2.  0.0.0  

2 2 0 0 0 

ENTER  (LINE  10)  NUMERICAL  SOlUTION  PARAMETERS  --  N(STD)  RDX  DT(STD)  TIWX  JMAX  --- 

20.  1.  1.  0.  0 -• 

20  1.  1.  0.  0 

ENTER  (LINE  11)  DIAGNOSTIC  PARAMETERS  --  IDIA6  IDCH  ID2  IPRZA  FLGPVE  — 

0.0.0. 0.0  

0 0 0 0 0. 

ENTER  (LINE  12)  PLOTTING  PARAMETERS  - IPLOTH  IPLOH  IPLOTH  IPLOTC  IPHAIT  --- 

0.0. 0.0.0  

0 0 0 0 0 

1 
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mmnmmmmmmmmmmtm  border  irrigation  flou  ttmtmtnntmmmntmmmmm 


EXAMPLE  tl 


ROUGHNESS  — 

HAWING  N (COEFFICIENT)  = 


INTAKE  CHARACTERISTICS  — 

K = 


A = 

B = 
C = 


08  JUN  83 


HYDRAULIC  PROPERTIES  OF 
DIMENSIONED 


CROP  AND  SOIL 
DIMENSIONLESS 


.0400  FT«(l/6)  = .0400  M»t(l/6)  EXPONENT  = 0. 

DIMENSIONLESS  DRAG  COEFFICIENT  Dt  = .lOOOEFOl 

IN  INFILTRATION  FORMU.A  Z = KtTBtA  4 BBT  f C 


1.0000  IN/HRBBA  = 2.5400  CH/HRttA 
.1291E400  IN/HINBBA  = .3279E401  HM/MINBBA 
.1389E-02  FT/SECIBA  = .4233E-03  M/SECBBA 


5000 

= .5000 

KB  = .1389E400 

A = .5000 

IN/HR 

= 0. 

CN/HR 

BB  = 0. 

IN 

= 0. 

CM 

CB  = 0. 

D I 

BORDER  G 
N E N S I 0 N E D 

E 0 H E T R Y 

DIMENSIONLESS 

18.21.01 
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FIELD  LENGTH  L =1200.0  FT  = 365.8  N LI  = .1410E+01 
AVG.  Fiat  aOPE  SOAVG  = .0005000  = .5000E-03  SOI  = .lOOOEFOl 
PLANE  SURFACE 

OPEN-END  BORDER 


IRRI6ATI0N-HANA6EKENT  PARAMETERS 

DIMENSIONED  DIMENSIONLESS 


REQUIRED  APPLICATION  = 1.50  IN  = 3.81  CM  ZREQI  = .2937Ei00 
UNIT  I)FLOU  RATE  = .2000  CFS/FOOT  =18.581  LITERS/SEC-M  QINI  = .lOOOETOl 
CUT-OFF  TINE  = 40.00  MIN  = .667  HR  TCOI  = .1325E401 


CHARACTERISTIC  PARAMETERS  IN  DIMENSIONLESS  FORMS 
CHMIACTERISTIC  DEPTH  BASED  ON  NORMAL  DEPTH  YO=YN>  XO=YO/SOt  T0=X0IY0/00 


DIMENSIONED  DIMENSIONLESS 

NORMAL  DEPTH  YN  = .4256  FT  = .1297  M = .1000E401 

FROUDE  NUMBER  AT  NORMAL  DEPTH  FN  = .1270E400 


CHARACTERISTIC  DEPTH  YO  = .4256E400  FT  = .1297E400  M 


CHARACTRISTIC  LENGTH 
CHAR'STIC  DISQMRGE 
CtMRACTERISTIC  TIME 


XO  = .8511E+03  FT 
00  = .2000E400  CFS/FOOT 
TO  = .3018E402  MIN 


= .2594E403  M 
= .1B58E+02  LITERS/SEC-M 
= .1811EF04  SEC 


TOP  = .2173E403  MIN  = .1304E405  SEC 

XOP  = .6128EI04  FT  = .1868EF04  N 

P = .7200E401 


NUMERICAL  SOLUTION  PARAMETERS 

ZERO-INERTIA  (EQUILIBRIUM)  SOLUTION 

DEPTH  MD  DISCHARGE  FUNCTIONS  ARE  TREATED  AS  FUUY  NONLINEAR 

MAXIMUM  ALLOWABLE  NUMBER  OF  ITERATIONS  JHAX  = 20 

WHEN  HATER  SURFACE  OF  COMPUTATIONAL  STREAM  IS  ESSENTIALLY  LEVa>  INFILTRATION  IS  ASSUMED  TO  BE  THAT  OF  LEVa  POOL(S) 
NUMBER  OF  CELLS  IN  STREM  = 20 

STANDARD  TIK  STEP  DT  = 1.00  MIN  = .017  W DTI  = .3313E-01 

DT=DT(STD)/N(STD)  FOR  FIRST  N(STD)  TINE  STEPS.  TIEN  OT=DT(STD> 

TNAX  = 780.00  MIN  = 13.0  HR  THAXI  = .2S84EF02 


ZERO-INERTIA  (EQUILIBRIUM)  NODa 
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POST-IRRIWTIW  SYNOPSIS 


K 

X 

TA 

1 

0. 

0. 

S 

.1410E+00 

.8171E-01 

9 

.2820EFOO 

.1897EW0 

13 

.4230EW0 

.3094E+00 

17 

.5640E+00 

.4349E+00 

21 

.7049E+00 

.5702EW0 

25 

.8459E+00 

.7080Ef00 

29 

.9869E+00 

.8495E+00 

33 

.1128E+01 

.9942EW0 

37 

.12«9Et01 

.lUlEFOl 

41 

.1410E+01 

.1410EI01 

.1291E+01 

TR 

Z 

.2887E+01 

.2360E400 

.3471E+01 

.2557EW0 

.4119E+01 

.2753E400 

.4584E401 

.2872EF00 

.4990E+01 

.29A4E400 

.5354E+01 

.3038EH)0 

.5A74EI01 

.3095E+00 

.^8E401 

.315SEI00 

.631SE401 

.3204E400 

.6603EF01 

.3248EF<)0 

.6761E401 

.3248EIOO 

0. 

INFILTRATION 


TA  MIN 

TR  MIN 

0. 

87.1 

2.5 

104.8 

5.7 

124.3 

9.3 

138.4 

13.2 

150.6 

17.2 

161.6 

21.4 

171.3 

25.6 

181.4 

30.0 

190.6 

34.5 

199.3 

39.0 

204.1 

PROFILE 


X METERS 

Z CM 

0. 

3.1 

36.6 

3.3 

73.2 

3.6 

109.7 

3.7 

146.3 

3.8 

182.9 

3.9 

219.5 

4.0 

256.0 

4.1 

292.6 

4.2 

329.2 

4.2 

365.8 

4.2 

3658EI03  0. 

X FEET 

Z IN 

0. 

1.21 

120.0 

1.31 

240.0 

1.41 

360.0 

1.47 

480.0 

1.51 

600.0 

1.55 

720.0 

1.58 

840.0 

1.61 

960.0 

1.64 

1080.0 

1.66 

1200.0 

1.66 

1200EF04  0. 

MEASURES  OF  MERIT  OF  COMPUTATION 


COMPUTATIONAL  VOLUME  BALANCE  VO=  .132SEf01  VZ=  .419&EH>0  VRO=  .9017ET00  RELATIVE  VOLUME  ERROR  = .2932E-02 

MAXIMUM  SIZE  OF  SAN  TEETH  IN  PROFILES  - ITM=  A RTM=  2 TOOTHM=  .138E+00  IRTM=228  KRTM=  3 RTOOTN=  .302E+01 


SYNOPSIS  OF  RESULTS  OF  IRRIGATION 


HOURS 

MINUTES 

DIMENSIONLESS 

T CO 

(TIIC  OF  CUT  OFF)  

40.00 

.1325E401 

T L 

(DURATION  OF  ADVMCE)  

38.97 

.1291E+01 

T R 

(TIME  RECESSION  STARTS  AT  UPSTREAM  END)  

87.15 

.2887EF01 

T FR 

(TIME  RECESSION  STARTS  AT  DOWNSTREAM  END)  

210.95 

.6989E+01 

T E 

(TINE  ALL  SURFACE  HATER  DISAPPEARS)  

210.95 

.6989EF01 

ENGLISH 

METRIC 

DIMENSIONLESS 

FTtt3/FT 

NU3/M 

V Q 

(APPLIED  VOLUME)  

44.594 

.1325EH)1 

V Z 

(INFILTRATED  VOLUS)  

14.119 

.4196EF00 

V RO 

(RUNOFF  VOLUME)  

30.345 

.9017EW0 

1 OF 

ENGLISH 

METRIC 

DIMENSIONLESS 

APPLIED 

DEPTH 

Y MAX  U 

(MAXIMUM  DEPTH  OF  SURFACE  STREAM  AHAINED  AT  UPPER  END)  

12.68  CM 

.9776E400 

Y MAX  D 

(MAXIMUM  DEPTH  OF  SURFACE  STREAM  ATTAINED  AT  LONER  END)  

0.  CM 

0. 

Y MAX 

(MAXIMUM  DEPTH  OF  SURFACE  STREAM  ATTAINED  AT  ANY  POINT)  

12.68  CM 

.9776E+00 

X YMAX 

(LOCATION  OF  POINT  OF  MAXIMUM  DEPTH)  

0.  N 

0. 

Z MIN 

(MINIMUM  DEPTH  OF  INFILTRATION)  

3.06  CM 

.2360EFOO 

25.1  Z 

Z MAX 

(MAXIMUM  DEPTH  OF  INFILTRATION)  

4.33  CM 

.3337EF00 

35.5  2 

Z LO 

(AVERAGE  LOH-OUARTER  DEPTH  OF  INFILTRATION)  

3.39  CM 

.2614EFOO 

27.8  Z 

Z REO 

(REQUIRED  DEPTH  OF  INFILTRATION)  

3.81  CM 

.2937E400 

31.3  Z 

Z Q 

(AVERAGE  (miED  DEPTH  VB/L)  

12.19  CM 

.9399EFOO 

100.0  Z 

Z AV6 

(AVERAGE  INFILTRATED  DEPTH  VZA)  

INCH 

3.86  CM 

.2976E+00 

31.7  2 

Z 

RO 
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(AVERAGE  DEPTH  OF  RUNOFF  VRO/L)  

3.27  INCH 

8.30  CM 

.6396E+00 

68 

.0  Z 

Z DP  ZREO 

(AVERAGE  DEPTH  OF  DEEP-PERCOLATION  VOLUME  BASED  ON  GIVEN  REQUIRED  DEPTH) 

.06  INCH 

.16  CM 

.1261E-01 

1.3  Z 

Z DP  ZHIN 

(AVERAGE  DEPTH  OF  DEEF-PERCOLATION  vaUME  BASED  ON  ZREQ  = ZMIN) 

.31  INCH 

.80  CM 

.6159E-01 

6.6  Z 

Z Df  ZLO 

(AVERAGE  DEPTH  OF  DEEF-PERCOLATION  VaUHE  BASED  ON  ZREO  = ZLO) 

.19  INCH 

.49  CM 

.3762E-01 

4.0  Z 

Z U ZREO 

(AVERAGE  OF  INFILTRATED  DEPTHS  LESS  THM  OR  EQUAL  TO  REQUIRED  DEPTH) 

1.46  INCH 

3.70  CM 

.28SOEIOO 

30.3  Z 

Z U ZMIN 

(AVERAGE  OF  INFILTRATED  DEPTHS  LESS  THAN  OR  EOU^  TO  ZMIN) 

1.21  INCH 

3.06  CM 

.2360EI00 

25.1  Z 

Z U ZLO 

(AVERAGE  OF  INFILTRATED  DEPTHS  LESS  THAN  OR  EQUAL  TO  ZLO) 

1.33  INCH 

3.37  CM 

.26OOEIO0 

27.7  Z 

UC  C 

(CHRISTIANSEN  UNIFORMITY  COEFFICIENT) 

.931 

UC  H 

(HSPA  UNIFORMITY  COEFFICIENT)  

.933 

DU 

(DISTRIBUTION  UNIFORMITY 

ZMIN/ZAVG)  

.793 

DU  LQ 

(LON-OUARTER  DISTRIBUTION  UNIFORMITY 

ZLQ/ZAVG)  

.879 

R? 

(RUNOFF  FRACTION  IN  PERCENT 

VRO/VQ  ) 

68.052 

BASED  ON 

BASED  ON 

BASED  ON 

GIVEN  ZREO 

ZREQ=ZNIN 

ZREO=ZLQ 

IE 

(IRRIGATION  EFFICIENCY 

ZU/ZO  ) 

25.11Z 

27.662 

UZ 

(USEFUL  FRACTION  OF  INFILTRATED  VOLUME 

ZU/ZAVG)  

.793 

.874 

SE 

(STORAGE  EFFICIENCY 

ZU/ZREQ)  

lOO.OOZ 

99.432 

AAP 

(PERCENT  OF  TOTAL  MEA  ADEQUATaY  IRRIGATED)  

lOO.OOZ 

88.082 

DR 

(DEFICIENCY  RATIO*  AVERAGE  DEFICIT  IN  UNDERIRRIGATED  MIEA* 

PERCENT  OF  ZREO)  8.11Z 

4.772 

CU-S  USED  FOR  THIS  RUN  = 39.044 


ttttiittmttmtmttttttttitBORDER-iRRiSATioN-FLoy  pROGRAHmtttmttttmtmmtimit 

08  JUN  83  18.24.02 

EU«>SED  CF  TINE  (SECONDS)  = 52.697 

CU-S  USED  = 158.196 

CU-S  REMAINING  = 541.804 


INTERACTIVE  USERS  --  POSITION  PAPER  FOR  PROMPTS  THEN  ENTER  (LINE  1)  1 FOR  A FRESH  START  — 

2 TO  CIMNGE  DESIGN  PARMETERS  - ZREOfOINiTCO  — 

3 TO  CHANGE  SOLUTION  MODES  - SOLMODfLIM(ODiDTNOD>ISUPZA>ZADMOD  — 

4 TO  CHM6E  SaUTIQN  PARAfSTERS  - NSTDfRDX>DTSTD>TMAXf  JHAX  — 

5 TO  CHANGE  LEVEL  OF  DIAGNOSTICS  AND/OR  PLOTTING  FLAGS  — 

0 TO  STOP  — 
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